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ABSTRACT 

We present analytic formulae that approximate the evolution of stars for a wide 
range of mass M and metallicity Z. Stellar luminosity, radius and core mass are given 
as a function of age, M and Z, for all phases from the zero-age main-sequence up 
to, and including, the remnant stages. For the most part we find continuous formulae 
accurate to within 5% of detailed models. These formulae are useful for purposes 
such as population synthesis that require very rapid but accurate evaluation of stellar 
properties, and in particular for use in combination with A^-body codes. We describe 
a mass loss prescription that can be used with these formulae and investigate the 
resulting stellar remnant distribution. 

Key words: methods: analytic - stars: evolution - stars: fundamental parameters - 
stars: Population II - galaxies: stellar content 



1 INTRODUCTION 

. The results of detailed stellar evolution calculations are re- 
' quired for applications in many areas of astrophysics. Ex- 
, amples include modelling the chemical evolution of galaxies, 
■ determining the ages of star clusters and simulating the out- 
comes of stellar collisions. As stellar evolution theory, and 
. our ability to model it, is continually being improved (the 
' treatment of convective overshooting and thermal pulses, 
for example) there is an ongoing need to update the results 
of these calculations. For a recent overview of problems in 
stellar evolution see Noels et al. (1995). 

As with all theories our understanding of stellar evo- 
lution must be tested against observations. One way to do 
this is to attempt to reproduce the findings of large-scale 
star surveys, such as the Bright Star Catalogue (Hoffleit 
1983) and the Hipparcos Catalogue (Ferryman et al. 1997), 
using population synthesis. The Hipparcos Catalogue is an 
excellent example of how improved observing techniques can 
initiate a re-evaluation of many aspects of stellar evolution 
theory (Baglin 1997; de Boer et al. 1997; Van Eck et al. 
1998). In order to make population synthesis statistically 
meaningful it is necessary to evolve a large sample of stars 
so as to overcome Foisson noise. If we synthesize n exam- 
ples of a particular type of star we have an error of ±^/n 
which means that for rarer stars often millions of possible 
progenitors are required to get a sufficently accurate sample. 
However, detailed evolution codes can take several hours to 



evolve a model of just one star. Thus it is desirable to gener- 
ate a large set of detailed models and present them in some 
convenient form in which it is relatively simple to utilise the 
results at a later stage. 

There are two alternative approaches to the problem 
of using the output of a series of stellar-evolution runs as 
data for projects that require them. One approach is to con- 
struct tables (necessarily rather large, especially if a range 
of metallicities and/or overshoot parameter is to be incor- 
porated) and interpolate within these tables. The other is to 
approximate the data by a number of interpolation formulae 
as functions of age, mass and metallicity. Both procedures 
have advantages and disadvantages (Eggleton 1996), so we 
have worked on both simultaneously. Stellar models have 
been available in tabular form for many years (Schaller et 
al. 1992; Charbonnel et al. 1993; Mowlavi et al. 1998). Steb 
lar populations cover a wide range of metallicity so the ideal 
is to have a set of models that cover the full range of pos- 
sible compositions and stellar masses. In a previous paper 
(Pols et al. 1998) we presented the results of stellar evolu- 
tion calculations for a wide range of mass and metallicity 
in tabular form. In the present paper we report on the re- 
sults of the second approach, construction of a set of single 
star evolution (SSE) formulae, thus expanding the work of 
Eggleton, Fitchett & Tout (1989) along the lines of Tout et 
al. (1996). It is more difficult in practice to find analytic ap- 
proximations of a conveniently simple nature for the highly 
non-uniform movement of a star in the Hertzsprung-Russell 
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diagram (HRD) than it is to interpolate in tables, but the 
resulting code is very much more compact and adaptable to 
the requirements of, for example, an A'^-body code (Aarseth 
1996) or variable mass loss. This is reinforced in the circum- 
stance where one wishes to include binary-star interactions, 
such as Roche-lobe overflow, common-envelope evolution, 
and magnetic braking with tidal friction, for example (Tout 
et al. 1997). 

In Section |^ we provide a brief overview of how stars 
behave as they evolve in time which introduces some of the 
terminology that we use and will hopefully facilitate the un- 
derstanding of this paper. Section H describes the detailed 
models from which the formulae are derived and justifies 
the inclusion of enhanced mixing processes. In Section ^ we 
outline the procedure to be used for generating the SSE 
package. The evolution formulae are presented in Section ^ 
for all nuclear burning phases from the main-sequence to the 
asymptotic giant branch. Our formulae are a vast improve- 
ment on the work of Eggleton, Fitchett & Tout (1989) not 
only due to the inclusion of metallicity as a free parameter 
but also because we have taken a great deal of effort to pro- 
vide a more detailed and accurate treatment of all phases of 
the evolution. Features such as main-sequence formulae that 
are continuous over the entire mass range and the modelling 
of second dredge-up and thermal pulses will be discussed. 
Section |^ discusses the behaviour of a star as the stellar en- 
velope becomes small in mass and outlines what happens 
when the nuclear evolution is terminated. We also provide 
formulae which model the subsequent remnant phases of 
evolution. In Section ^ we describe a comprehensive mass 
loss algorithm which can be used in conjunction with the 
evolution formulae, as well as a method for modelling stellar 
rotation. Various uses for the formulae and future improve- 
ments are discussed in Section ^ along with details of how 
to obtain the formulae in convenient subroutine form. 



2 STELLAR EVOLUTION OVERVIEW 

A fundamental tool in understanding stellar evolution is the 
Hertzsprung-Russell diagram (HRD) which provides a cor- 
relation between the observable stellar properties of lumi- 
nosity, L, and effective surface temperature, Teff. Figure ^ 
shows the evolution of a selection of stars in the HRD from 
the zero-age main-sequence (ZAMS), where a star adjusts 
itself to nuclear burning equilibrium, until the end of their 
nuclear burning lifetimes. As stars take a relatively short 
time to reach the ZAMS all ages are measured from this 
point. The length of a stars life, its path on the HRD and 
its ultimate fate depend critically on its mass. 

Stars spend most of their time on or near the main- 
sequence (MS) burning hydrogen to produce helium in their 
cores. To first order, the behaviour of a star on the MS can 
be linked to whether it has a radiative or convective core. 
Stars with M ^ 1.1 Mq have radiative cores while in higher 
mass stars a convective core develops as a result of the steep 
temperature gradient in the interior. During core hydrogen 
burning on the MS, low-mass stars will move upwards in L 
and to higher Tctt on the HRD while higher mass stars will 
also move upwards in L but to a region of lower Tgff . The MS 
evolution will end when the star has exhausted its supply of 
hydrogen in the core. Low-mass stars will continue expand- 



ing as they evolve off the MS but for higher mass stars with 
convective cores the transition is not so smooth. Owing to 
mixing in the core there is a sudden depletion of fuel over 
a large region which leads to a rapid contraction over the 
inner region at core hydrogen-exhaustion. This causes the 
hydrogen-exhausted phase gap, or MS hook, which occurs 
on a thermal timescale. The different features of MS evo- 
lution are illustrated by comparing the evolution tracks for 
the 1.0 Mq and 1.6 Mq stars in Figure |l} 

The immediate post-MS evolution towards the right in 
the HRD occurs at nearly constant luminosity and is very 
rapid. For this reason very few stars are seen in this phase, 
and this region of the HRD is called the Hertzsprung gap 
(HG), or the sub-giant branch. During this HG phase the ra- 
dius of the star increases greatly causing a decrease in T^s- 
For cool envelope temperatures the opacity increases causing 
a convective envelope to develop. As the convective envelope 
grows in extent the star will reach the giant branch (GB) 
which is the nearly vertical line corresponding to a fully 
convective star, also known as the Hayashi track. All stars 
ascend the GB with the hydrogen-exhausted core contract- 
ing slowly in radius and heating while the hydrogen-burning 
shell is eating its way outwards in mass and leaving behind 
helium to add to the growing core. As the stars move up 
the GB convection extends over an increasing portion of the 
star. The convective envelope may even reach into the previ- 
ously burnt (or processed) regions so that burning products 
are mixed to the surface in a process called dredge-up. 

Eventually a point is reached on the GB where the core 
temperature is high enough for stars to ignite their central 
helium supply. For massive stars, M ^ 2.0 Mq , this takes 
place gently. When core helium burning (CHeB) begins the 
star descends along the GB until contraction moves the star 
away from the fully convective region of the HRD and back 
towards the MS in what is called a blue loop. During CHeB, 
carbon and oxygen are produced in the core. Eventually core 
helium is exhausted and the star moves back to the right in 
the HRD. The size of the blue loop generally increase with 
mass, as can be seen by comparing the 4.0 Mq and 10.0 Mq 
tracks in Figure Lower mass stars have degenerate helium 
cores on the GB leading to an abrupt core-helium flash at 
helium ignition (Hel). The star then moves down to the 
zero-age horizontal branch (ZAHB) very quickly. The initial 
position of a star along the ZAHB depends on the mass of 
the hydrogen-exhausted core at the time of ignition and on 
the mass in the overlying envelope. Those stars with lower 
mass, ie. shallower envelopes, appear bluer because there is 
less mass to shield the hot hydrogen burning shell. It is also 
possible for stars of very high mass, Af J; 12.0 Mq, to reach 
high enough central temperatures on the HG for helium to 
ignite before reaching the GB. The 16.0 Mq star in Figure ^ 
is such an example. As a result these stars by-pass the GB 
phase of evolution. 

Evolution after the exhaustion of core-helium is very 
similar to evolution after core-hydrogen exhaustion at the 
end of the MS. The convective envelope deepens again so 
that the star once more moves across towards the Hayashi 
track to begin what is called the asymptotic giant branch 
(AGB). On the AGB the star consists of a dense core com- 
posed of carbon and oxygen surrounded by a helium burn- 
ing shell which adds carbon to the degenerate core. Initially 
the H-burning shell is extinguished so that the luminosity is 
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supplied exclusively by the He-burning shell; characterizing 
the early AGB (EAGB) phase. If the star is massive enough 
the convective envelope can reach into the H-exhausted re- 
gion again (second dredge-up). When the He- burning shell 
catches up with the H-rich envelope the H-shell reignites 
and the two grow together with the H-burning shell sup- 
plying most of the luminosity. During the following phase 
the helium shell is unstable, which can cause a helium shell 
flash in which the helium shell will suddenly release a large 
amount of luminosity. The energy released in the flash ex- 
pands the star resulting in the hydrogen shell cooling so 
much that it is extinguished. Convection once again reaches 
downward past the dead hydrogen shell. This mixes helium 
to the surface, as well as carbon that was mixed out of the 
helium shell by flash-driven convection. As the star subse- 
quently contracts the convection recedes and the hydrogen 
shell re-ignites but has now moved inwards in mass due to 
the envelope convection. This process is called third dredge- 
up. The star continues its evolution up the AGB with the 
hydrogen shell producing almost all of the luminosity. The 
helium shell flash can repeat itself many times and the cycle 
is known as a thermal pulse. This is the thermally pulsing 
asymptotic giant branch (TPAGB). 

The stellar radius can grow to very large values on the 
AGB which lowers the surface gravity of the star, so that the 
surface material is less tightly bound. Thus mass loss from 
the stellar surface can become significant with the rate of 
mass loss actually accelerating with time during continued 
evolution up the AGB. Unfortunately, our understanding of 
the mechanisms that cause this mass loss is poor with pos- 
sible suggestions linking it to the helium shell flashes or to 
periodic envelope pulsations. Whatever the cause, the influ- 
ence on the evolution of AGB stars is significant. Mass loss 
will eventually remove all of the stars envelope so that the 
hydrogen burning shell shines through. The star then leaves 
the AGB and evolves to hotter T^s at nearly constant lumi- 
nosity. As the photosphere gets hotter the energetic photons 
become absorbed by the material which was thrown off while 
on the AGB. This causes the material to radiate and the star 
may be seen as a planetary nebula. The core of the star then 
begins to fade as the nuclear burning ceases. The star is now 
a white dwarf (WD) and cools slowly at high temperature 
but low luminosity. 

If the mass of the star is large enough, M <: 7M0, the 
carbon-oxygen core is not degenerate and will ignite car- 
bon as it contracts, followed by a succession of nuclear re- 
action sequences which very quickly produce an inner iron 
core. Any further reactions are endothermic and cannot con- 
tribute to the luminosity of the star. Photodisintegration of 
iron, combined with electron capture by protons and heavy 
nuclei, then removes most of the electron degeneracy pres- 
sure that was supporting the core and it begins to collapse 
rapidly. When the density becomes large enough the inner 
core rebounds sending a Shockwave outwards through the 
outer layers of the star that have remained suspended above 
the collapsing core. As a result the envelope of the star is 
ejected in a supernova (SN) explosion so that the AGB is 
effectively truncated at the start of carbon burning and the 
star has no TPAGB phase. The remnant in the inner core 
will stablise to form a neutron star (NS) supported by neu- 
tron degeneracy pressure unless the initial stellar mass is 



large enough that complete collapse to a black hole (BH) 
occurs. 

Stars with AI <: 15 M0 are severely affected by mass 
loss during their entire evolution and may lose their en- 
velopes during CHeB, or even on the HG, exposing nuclear 
processed material. If this occurs then a naked helium star 
is produced and such stars, or stars about to become naked 
helium stars, may be Wolf-Rayet stars. Wolf-Rayet stars are 
massive objects which are found near the MS, are losing 
mass at very high rates and show weak, or no, hydrogen 
lines in their spectra. Luminous Blue Variables (LBVs) are 
extremely massive post-MS objects with enormous mass loss 
rates in a stage of evolution just prior to becoming a Wolf- 
Rayet star. Naked helium stars can also be produced from 
less massive stars in binaries as a consequence of mass trans- 
fer. 

Variations in composition can also affect the stellar evo- 
lution timescales as well as the appearance of the evolution 
on the HRD, and even the ultimate fate of the star. A more 
detailed discussion of the various phases of evolution can be 
found throughout this paper. 



3 STELLAR MODELS 

The fitting formulae are based on the stellar models com- 
puted by Pols et al. (1998). They computed a grid of evo- 
lution tracks for masses M between 0.5 and 50 Mq and 
for seven values of metallicity, Z = 0.0001, 0.0003, 0.001, 
0.004, 0.01, 0.02 and 0.03. They also considered the prob- 
lem of enhanced mixing such as overshooting beyond the 
classical boundary of convective instability. Its effect was 
modelled with a prescription based on a modification of the 
Schwarzschild stability criterion, introducing a free param- 
eter 5ov (which differs from the more commonly used pa- 
rameter relating the overshooting distance to the pressure 
scale height; see Pols et al. 1998 for details). The tracks com- 
puted with a moderate amount of enhanced mixing (given 
by 5ov = 0.12 and labeled the OVS tracks by Pols et al. 
1998) were found to best reproduce observations in a series 
of sensitive tests involving open clusters and ecliping bina- 
ries (see Schroder, Pols & Eggleton 1997; Pols et al. 1997, 
1998). We consequently use these OVS tracks as the data to 
which we fit our formulae. 

For each Z, 25 tracks were computed spaced by approx- 
imately 0.1 in logM, except between 0.8 and 2.0 Mq where 
four extra models were added to resolve the shape of the 
main-sequence which changes rapidly in this mass range. 
Hence we dispose of a database of 175 evolution tracks, each 
containing several thousand individual models. 

A subset of the resulting OVS tracks in the HRD are 
shown in Fig. for Z = 0.02 and Fig | for Z = 0.001. 
The considerable variation of model behaviour introduced 
by changes in metallicity is illustrated by Fig. ^. Detailed 
models of the same mass, M = 6.35 Mq, are shown on the 
HRD for three different metallicities, Z = 0.0001, 0.001 and 
0.02. Not only does a change in composition move the track 
to a different position in the HRD but it also changes the 
appearance of each track, as can be seen by considering the 
extent of the hook feature towards the end of the main se- 
quence and the blue loops during core helium burning. Fur- 
thermore, the Z — 0.0001 model ignites helium in its core 
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Figure 1. Selected OVS evolution tracks for Z = 0.02, for masses 
0.64, 1.0, 1.6, 2.5, 4.0, 6.35, 10, 16, 25 and 40 Mq. 
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Figure 2. Same as Fig. |l| for Z = 0.001. The 1.0 post He 
flash track has been omitted for clarity. 



while on the Hertzsprung gap as opposed to the other mod- 
els which evolve up the giant branch before reaching a high 
enough core temperature to start helium burning. In addi- 
tion the nuclear burning lifetime of a star can change by as 
much as a factor of 2 owing to differences in composition, 
as shown in Fig. ^ for a set of 2.5 Mq models. This em- 
phasizes the need to present the results of stellar evolution 
calculations for an extensive range of metallicity. 

Mass loss from stellar winds was neglected in the de- 
tailed stellar models, mainly because the mass loss rates are 
uncertain by at least a factor of three. We do include mass 
loss in our analytic formulae in an elegant way, as will be de- 
scribed in Section which allows us to experiment easily 
with different mass loss rates and prescriptions. 



4 PROCEDURE 

We assign each evolution phase an integer type, k, where: 

— MS star M 0.7 deeply or fully convective 

1 = MS star M ^ 0.7 

2 — Hertzsprung Gap (HG) 

3 = First Giant Branch (GB) 

4 = Core Helium Burning (CHeB) 

5 — Early Asymptotic Giant Branch (EAGB) 



6 = Thermally Pulsing Asymptotic Giant Branch (TPAGB) 

7 = Naked Helium Star MS (HeMS) 

8 — Naked Helium Star Hertzsprung Gap (HeHG) 

9 = Naked Hehum Star Giant Branch (HeCB) 

10 = Helium White Dwarf (He WD) 

11 = Carbon/Oxygen White Dwarf (CO WD) 

12 = Oxygen/Neon White Dwarf (ONe WD) 

13 = Neutron Star (NS) 

14 = Black Hole (BH) 

15 — massless remnant, 

and we divide the MS into two phases to distinguish be- 
tween deeply or fully convective low-mass stars and stars of 
higher mass with little or no convective envelope as these 
will respond differently to mass loss. 

To begin with we take different features of the evolution 
in turn, e.g. MS lifetime, ZAHB luminosity, and first try to 
fit them as / (M) for a particular Z in order to get an idea of 
the functional form. We then extend the function to g (M, Z) 
using / (Af ) as a starting point. In this way we fit formulae 
to the end-points of the various evolutionary phases as well 
as to the timescales. We then fit the behaviour within each 
phase as h (t, M, Z), e.g. Lms {t, M, Z). 

As a starting point we take the work of Tout et 
al. (1996) who fitted the zero-age main-sequence luminos- 
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Figure 3. Detailed OVS evolution tracks for M = 6.35 Mq, for 
metallicities 0.0001, 0.001 and 0.02. 
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Figure 4. Radius evolution as a function of stellar age for 
M = 2.5 Mq, for metallicities 0.0001, 0.001 and 0.02. Tracks are 
from the detailed models and run from the ZAMS to the point of 
termination on the AGB. 



ity (LzAMs) and radius (i?zAMs) as a function of M and Z. 
Their aim, as is ours, was to find simple computationally 
efficent functions which are accurate, continuous and differ- 
entiable in M and Z, such as rational polynomials. This is 
acheived using least-squares fitting to the data after choos- 
ing the initial functional form. In most cases we determine 
the type of function, the value of the powers, and the num- 
ber of coefficients to be used, simply by inspecting the shape 
of the data, however in some cases, such as the luminosity- 
core-mass relation on the giant branch, the choice will be 
dictated by an underlying physical process. For the ZAMS, 
accuracy is very important because it fixes the star's posi- 
tion in the HRD. Tout et al. (1996) acheived Lzams accurate 
to 3% and i?zAMS accurate to 1.2% over the entire range. 
For the remainder of the functions we aim for RMS errors 
less than 5% and preferably a maximum individual error 
less than 5% although this has to be relaxed for some later 
stages of the evolution where the behaviour varies greatly 
with Z but also where the model points are more uncertain 
owing to shortcomings in stellar evolution theory. 



5 FITTING FORMULAE 

In this section we present our formulae describing the evo- 
lution as a function of mass M and age t. The explicit Z- 
dependence is in most cases not given here because it would 
clutter up the presentation. This ^-dependence is implicit 
whenever a coefficient of the form a„ or b„ appears in any of 
the formulae. The explicit dependence of these coefficients 
on Z is given in the Appendix. Coefficients of the form Cn, 
whose numerical values are given in this section, do not de- 
pend on Z. 

We adopt the following unit conventions: numerical val- 
ues of mass, luminosity and radius are in solar units, and 
values of timescales and ages are in units of lO'' yr, unless 
otherwise specified. 

We begin by giving formulae for the most important 
critical masses, Mhook (the initial mass above which a hook 
appears in the main-sequence), Mhgp (the maximum initial 
mass for which He ignites degenerately in a helium fiash) 
and A/pGB (the maximum initial mass for which He ignites 
on the first giant branch). Values for these masses are given 
in Table 1 of Pols et al. (1998) estimated from the detailed 
models for 7 metallicities. These values can be accurately 
fitted as a function of Z by the following formulae, where 
C = log(Z/0.02): 



Mhook = 1.0185 -I- 0.16015C + 0.0892C , 
A/hof = 1.995 + 0.25C + 0.087C^ 



AfpGB = 



13.048 (z/omy 



(1) 

(2) 
(3) 



1-H 0.0012 (0.02/Z)^'^^' 

Based on the last two critical masses, we make a distinction 
into three mass intervals, which will be useful in the later 
descriptions: 

(i) low-mass (LM) stars, with M < MhbF, develope de- 
generate He cores on the GB and ignite He in a degenerate 
fiash at the top of the GB; 

(ii) intermediate-mass (IM) stars, with MhbF < M < 
Mfgb , which evolve to the GB without developing degener- 
ate He cores, also igniting He at the top of the GB; 



© 1999 RAS, MNRAS 000, 



6 J. R. Hurley, O. R. Pols and C. A. Tout 



(iii) high-mass (HM) stars, with M > Mfgb , ignite He in 
the HG before the GB is reached, and consequently do not 
have a GB phase. 

Note that this definition of IM and HM stars is different from 
the more often used one, based on whether or not carbon 
ignites non-degenerately. 

5.1 Main-sequence and Hertzsprung gap 

To determine the base of the giant branch (BGB) we find 
where the mass of the convective envelope Mce first exceeds 
a set fraction of the envelope mass Me as Mce increases on 
the HG. From inspection the following fractions 

2 

— 1 
5 

1 . 
3^ 

generally give a BGB point corresponding to the local mini- 
mum in luminosity at the start of the GB. We define helium 
ignition as the point where Lhc = O.OIL for the first time. 
For HM stars this will occur before the BGB point is found, 
ie. no GB, and thus we set tsGB = teei for the sake of defin- 
ing an end-point to the HG, so that BGB is more correctly 
the end of the HG (EHG) as this is true over the entire mass 
range. 

The resultant lifetimes to the BGB are fitted as a func- 
tion of M and Z by 

ai + asM* + asAf^'^' + M'^ 



Mce 



;Me M <, MneF 
-Me M ^ MhoF 



tsGB = 



a^M^ + asM'^ 



(4) 



Figure g shows how eq. (j^ fits the detailed model points for 
Z = 0.0001 and 0.03 which are the metallicities which lead 
to the largest errors. Over the entire metallicity range the 
function gives a rms error of 1.9% and a maximum error of 
4.8%. In order that the time spent on the HG will always be 
a small fraction of the time taken to reach the BGB, even 
for low-mass stars which don't have a well defined HG, the 
MS lifetimes are taken to be 

tms = niax (thook, xtBCs) , (5) 
where thook = M*bgb and 

X = max (0.95, min (0.95 - 0.03 (C + 0.30103) , 0.99)) (6) 



jj, — max ( 0.5, 1.0 — 0.01 max 



(7) 



Note that jj, is ineffective for M < Afhook, ie. stars without 
a hook feature, and in this case the functions ensure that 

X > jj.. 

So we now have defined the time at the end of the MS, 
tus and the time taken to reach the start of the GB (or end 
of the HG), tBGB such that 



t : 0.0 ^ tMS 



MS evolution 
HG evolution. 



The starting values for L and R are the ZAMS points 
fitted by Tout et al. (1996). We fit the values at the end of 
the MS, Ltms and -Rtms, as well as at the end of the HG, 



Lehg ~ 



Lbcb M < Mfgb 
Luoi M > Mfgb 



Z = 0.0001 
Z = 0,03 



log(M/Mg) 

Figure 5. Time taken to reach the base of the giant branch as 
a function of stellar mass as given by eq. (^), shown against the 
detailed model points, for Z = 0.0001 and 0.03 which give the 
worst fit of all the metallicities. The maximum error over the 
entire metallicity range is 4.8% and the RMS error is 1.9%. 



Rehg 



Rgb {Lbgb) M < Mfgb 



i?Hoi M > Mfgb 
The luminosity at the end of the MS is approximated 



by 



Ltms = 



aiiM^' + a^^M^ + ai3M"i«+i'* 



(8) 



ai4 + aif,M'^ + M'Jie 

with ttie ~ 7.2. This proved fairly straightforward to fit but 
the behaviour of -Rtms is not so smooth and thus requires 
a more complicated function in order to fit it continuously. 
The resulting fit is 

ais + aiaM°-^^ 



Rtms 



R 



CiM^ + assM^^e + a24M°-^o+^'^ 



(9) 



M> M,, (9a) 



a25 + MS 

with straight-line interpolation to connect eqs. (^) and (|^) 
between the endpoints, where 

M, = an + 0.1 , an ^ 1.4 

and ci = -8.672073 x 10~^ Oai « 1.47, a22 ~ 3.07, ase ~ 
5.50. Note that for low masses, M < 0.5, where the function 
is being extrapolated we add the condition 

-Rtms = max (-Rtms, 1.5-Rzams) 

to avoid possible trouble in the distant future. 

The luminosity at the base of the GB is approximated 

by 

a29 + a^oMCs + Ma32 ' ^^^> 

with C2 = 9.301992, C3 = 4.637345, O31 4.60 and aga ~ 
6.68. The description of -LhoI, -Rgb and 7?hcI is given in 
later sections. 



5.1.1 Main-sequence evolution 

On the MS we define a fractional timescale 
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Figure 6. Luminosity evolution on the main-sequence for a typ- 
ical detailed model with a hook feature (star points) decomposed 
into two functions: a smooth polynomial (solid line) and a pertu- 
bation function (dash-dot line). 



Lb{t) for a typical detailed model is illustrated by Fig. ^ 
The luminosity pertubation is approximated by 



r 0.0 



AL : 



B 



M - Mhook 



133 — Afhook- 

/ 034 ase \ 

AM"35'Ma37j 



M < Mhook 

Mhook < M < 033 (16) 
M > 033 



where B — AL (033), 1.25 < 033 < 1.4, 035 ~ 0.4 and 037 
0.6. 

The radius pertubation is approximated by 



( 0.0 



043 



/ M- Mhook y 

V 042 — Mhook / 



AT? = < 



043 + {B — 043) 

038 + a39M^-^ 
I a4oM3 -I- M"4i 



M - O42 



2.0 - 042 
- 1.0 



0.44 



M < Mhook 
Mhook < M < 04 
042 < M < 2.0 
M > 2.0 



(17) 



t 



(11) 



As a star evolves across the MS its evolution accelerates 
so that it's possible to model the time dependence of the 
logarithms of the luminosity and radius by polynomials in 
r. Luminosity is given by 



log 



izAMS 



/J^r" + (log ■ 



jtms 



- Q-L 



-AL(rl-rl 



(12) 



and radius by 



log 



Rvis (t) 

-RzAMS 

+ (log 



= Qht -I- Prt +-yT + 

RtMS r, \ 3 A D / 3 

— an — Pn —'y]T — AR ( ri 



(13) 



Rz 



AMS 



where 

n = min (1.0, t/ihook) 
T2 — max I 0.0, min ( 1.0 



t — (1.0 — e) fhook 

t ^hook 



(14) 
(15) 



for e = 0.01. 

We add AL and AR as pertubations to the smooth 
polynomial evolution of L and R in order to mimic the hook 
behaviour for M > A/hook- In effect we have 

Lms it) = La (t) /Lb (t) 

where La (t) is a smooth function describing the long-term 
behaviour of Lms (i) and Lt (t) is another smooth function 
describing short-term pertubations where 

logLt, {t) = AL (n -r|) 

and the action of T2 acheives a smooth transition over 
At — ethook- This decomposition of L{t) into La{t) and 



where B = AR{M = 2.0), 041 « 3.57, 1.1 < O42 < 1.25 and 
044 ~ 1.0. 

The exponent 77 = 10 in eq. (|l|) unless Z < 0.0009 
when it is given by 



10 M < 1.0 
20 M > 1.1 



(18) 



with linear interpolation between the mass limits. 

The remaining functions for this section are those that 
describe the behaviour of the coefBcients in eqs. (^2|) and 
(p^. The fact that these can appear messy and complicated 
in places reflects rapid changes in the shape of the L and 
R evolution for the detailed models as a function of Al as 
well as Z. This is illustrated in Figs. ^, ^, |^ and |l^ which 
also show the tracks derived from these functions, exhibit- 
ing that our efforts have not been in vain. The fitting of the 
coefficients is also complicated by the sensitivity of eqs. ( |l2| ) 
and (^) to small changes in the values of the coefflcients. 
Ideally we would like all the functions to be smooth and 
different iable across the entire parameter space but in some 
places this has to be sacrificed to ensure that the position of 
all the fitted tracks on the HRD is as accurate as possible. 
This is deemed necessary as the main use of the functions 
is envisaged to be the simulation of Colour-Magnitude dia- 
grams for comparison with observations. 

The luminosity a coefficient is approximated by 



O48 



O45 -I- a4BM 
MO-4 a47Mi-' 



M > 2.0, 



(19) 



where 04s ~ 1.56, and then 
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Figure 7. Luminosity evolution on the main-sequence as given 
by eq. (solid line) and from the detailed models (points) for 
selected masses with a metallicity of 0.001. 



M < 0.5 

6149 + 5.0 (0.3 - aio) (M - 0.5) 

0.5 < M < 0.7 

0.3 + (aso - 0.3) {M - 0.7) / - 0.7) 
0.7 < M < a52 

^50 + (flsi — tlso) {M ~ 052) / (Cl53 — 052) 

052 < M < ass 

asi + (B - asi) (M - 053) / (2.0 - 053) 
flBs < M < 2.0 



where B = (M ^ 2.0). 

The luminosity /3 coefficient is approximated by 

(Bl = max (0.0, a54 - a^sM""^'') 

where Ose ~ 0.96. Then if M > 057 and Pl > 0.0 

/3l = max (0.0, B - 10.0 (A/ - a^r) B) 

where B = Pl {M = 057) and 1.25 < 057 < 1.4. 
The radius a coefficient is approximated by 



(Ma) 



066 < A/ < fle? 



where Oee ~ 1.4 and a67 « 5.2, and then 



(21) 




Figure 8. Same as Fig. ^ for Z = 0. 



02. 



( 062 



M < 0.5 



062 + (063 - O62) (M - 0.5) /0.15 

0.5 <M< 0.65 

063 + (o64 - 063) (M - 0.65) / (oes - 0.65) 

0.65 < M < 06S 

064 + {B — O64) (M — Oes) / (O66 — 068) 

Oes < M < 066 



. C + 065 {M - O67) 



067 < M 



where B = an {M — Ogg) and C = a_R. (Af = 0^7). 

The radius /3 coefficient is approximated by Pr — — 1, 
where 



P'r 



aeaM' 



3.5 



O70 + Af"71 

with 071 ~ 3.45, and then 
( 1.06 



2.0 < M < 16.0 , 



(22) 



(20) p'^^ ^ 



M < 1.0 



(P* 



1.06 + (072 - 1.06) {M - 1.0) / (074 - 1.06) 

1.0 < M < 074 

072 + (B - 072) {M - 074) / (2.0 - 074) 

074 < M < 2.0 

C + 073 (A4" - 16.0) 16.0 < M 

where B = P'r (M = 2.0), C = (Af = 16.0). 

If M > + 0.1 then 7 = 0.0 where O75 ~ 1.25. Other- 
wise 

076 + O77 {M — 073)"^^ 

M - 1.0 °- 



7 ' 



B + (080 - B) 



.ars - 1.0 
C - 10.0 (A/ - 075) C 
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M = 0,80 
Z = 0.001 




M = 1.00 
Z = 0.020 




M = 1.60 
Z = 0.001 




Figure 9. Radius evolution on the main-sequence as given by 
eq. ( p^ ) (solid line) and from the detailed models (points) for 
selected masses with a metallicity of 0.001. 



where B = -y {M = 1.0) and C = ago unless 075 < 1.0 when 
C = B. Note we must always double-check that 7 > 0.0. 

Following Tout et al. (1997) we note that low-mass MS 
stars can be substantially degenerate below about 0.1 Mq so 
we take 



Rms = max (Ems, 0.0258 (1-0 + Xf'^ hr^'^) 
for such stars. 

5.1.2 Herizsprung gap evolution 

During the HG we define 

t — tms 
T — . 

tBGB — ^MS 

Then for the luminosity and radius we simply take 

/ I/EHG \ ^ 



(24) 



Lug ~ ^TMS 
Rug ~ Rtms 



\ I/TMS ' 

/ RmiG \ ^ 

V i?TMS / 



(25) 



(26) 
(27) 



On the MS we don't consider the core to be dense 
enough with respect to the envelope to actually define a 
core mass, ie. Mc.ms = 0.0. The core mass at the end of the 
HG is 

r Mc,GB {L = Lbgb) M < Mhcf 

Mc,EHG = < Mc,BGB MhoF < M < MpGB (28) 

I Mc.Hci M > Mfgb , 

where Mc,gb, Afc bgb and Afc,Hci will be defined in Sec- 
and At the beginning of the HG we set 
/oMcEHG, where 

^ 1.586 + MS-^s 
^ 2.434 + 1.02M5-25 ' ^ ^> 

and simply allow the core mass to grow linearly with time 
so that 



tions 5.2 

Mc,TMS 




Figure 10. Same as Fig. |9| for Z = 0. 



Mc.HG = [{l-~r)p + t] Mc.ehg • 



02. 



(30) 



If the HG star is losing mass (as described in Section 7.1) it 
is necessary to take Mc,hg as the maximum of the core mass 
at the previous timestep and the value given by eq. 

5.2 First giant branch 

The evolution along the first giant branch (GB) can 
elled, foUowing Eggleton, Fitchett & Tout (1989), 
power-law core mass-luminosity relation, 

The evolution is then determined by the growth of 
mass as a result of H burning which, in a state of 
equilibrium, is given by 

L = EXeMc A/c = AnL 

where 

Xe — envelope mass fraction of hydrogen, 

E — the specific energy release and 
Ah ~ hydrogen rate constant. 

Thus 



be mod- 
using a 

(31) 

the core 
thermal 

(32) 



,^ -AuDM^ 
at 

which upon integration gives 
= [{p - 1) AnD {ti^f - t)]rh 

or 

L = D[{p~l)AnD{ti^t-t)]'^ 



(33) 

(34) 
(35) 



so that the time evolution of either Mc or L is given and 
we can then simply find the other from the M^-L relation. 
Also, when L = Lbgb we have t — tsGB which defines the 
integration constant 
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tini ~ i 



1 / D y 

^^''^ + AuD (p - 1) [l^J 



(36) 



Now as noted in Tout et al. (1997), the single power- 
law L oc Mc is a good approximation to the evolution for 
small Mc but the relation flattens out as Mc approaches the 
Chandrasekhar mass A'Ich- They expanded the relation to 
consist of two power-law parts. We use an improved form 
which, albeit somewhat more ad hoc, follows much better 
the actual time evolution along the GB. Our Mc-L relation 
has the form 



L = min {BM,", DM,") {q < p), 



(37) 



so that the first part describes the high-luminosity end and 
the second the low-L end of the relation with the two cross- 
ing at 



(i) 



(38) 



The parameters B, D, p and q are constants in time for 
each model and indeed are constant in mass for Af < MhgF- 
For M > MhbF it is necessary to introduce a dependence 
on initial mass so that we actually have a M^-L-M relation. 
The only region in the Mc-L parameter space where we find 
that a Z-dependence is required is in the value of D for 
M < AiucF- The parameters are 



p = 



6 M < Mhof 

5 M> 2.5 

3 M < Mhof 

2 M > 2.5 



B 



logD 



(3 X lO", 500 + 1.75 X 10* M° '') 

r 5.37 + 0.135C [= Do] M < Mhof 

max (-1.0, 0.975 A) - 0.18M, 0.5i3o - 0.06M) 
1 M> 2.5 



with linear interpolation over the transition region, MhbF < 
M < 2.5, in order to keep the parameters continuous in M. 
Thus isochrones constructed with these functions will not 



give a discontinuity on the GB. The behaviour of eq. (|37 
is shown in Fig. [Ill as the fit to selected model points (note 
how the relation flattens out as the luminosity increases). 
Equation d34h now becomes 



Mc,GB = 



[{p-l)AHD(tinl,l-t)]— t<t^ 



[(g - 1) AhB (tinf,2 

for Ibgb <t< tuei, where 

^ {p - 1) AhD (lbgb) 

(tinf.l — ^BGb) ( ) 



t > to: 



tini, 



tsGB 
iinf.l 



p-1 

D \ — 



tini,2 ~ 



1 / B\ 

+ (g-i)AHS virJ 



B \ V 



(39) 

(40) 
(41) 
(42) 



The G B en ds ai t — tuci, corresponding to L = LhoI (see 
Section 5^), given by 




M = 1.6 Z = 0.0003 
M = 2.5 Z = 0.0200 
M = 1.0 Z = 0.0200 



1 



1.2 



Figure 11. Relation between core mass and luminosity on the 
giant branch showing the fit to points taken from selected detailed 
models given by eq. (solid lines). 



M 


log^ii 


tHoI-*BGB 
tBGB 


tBGB/Myr 


1.0 
2.0 
5.0 


-4.8 
-4.1 
-3.4 


6.4 X 10-2 
2.0 X 10-2 
2.4 X 10-3 


10* 
10^ 

102 



Table 1 . A selection of values for the mass dependent hydrogen 
rate constant with approximate timescales also listed. 



iiiif,i 



tu 



p-1 

1 ( B \ " , 



(43) 



The value used for An depends on whether we take 
the PP chain or the CNO cycle as the hydrogen burning 
mechanism with the CNO cycle being the most likely on the 
GB. Now 

E = ecNo/m (i?e*) ^ 6.018 x 10^** erg/g 
thus 

Au = {EX^y^ = 2.37383 x 10"^^ g/erg 

=> Au^ 1.44 X lO"-"' Mq LQ~^Myr~^ , 

i.e. log Ah = —4.84. In practice there are small deviations 
from thermal equilibrium which increase with stellar mass. 
As the value of Au fixes the rate of evolution on the GB and 
thus the GB timescale it is important for it to be accurate 
especially if we want to use the formulae for population syn- 
thesis. We find that the detailed models are best represented 
if we introduce a mass dependant An, ie. A'^, where 

log A'u = max (-4.8, min (-5.7 + 0.8M, -4.1 + 0.14Af)) . 

Some representative values of A'^ as a function of initial 
stellar mass are shown in Table 1 along with approximate 
values for the GB lifetime and the time taken to reach the 
GB. 

Evolution on the GB actually falls into two fairly dis- 
tinct categories depending on whether the initial mass of the 
star is greater than or less than MhcF- If M < A4"hoF then 
the star has a degenerate helium core on the GB which grows 
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according to the Mc,gb relation derived from eq. @. When 
helium ignites at the tip of the GB it does so degenerately 
resulting in the helium flash. However, for IM stars on the 
GB, M > MhgF, the helium core is non-degenerate and the 
relative time spent on the GB is much shorter and thus the 
models show that Mc,gb is approximately constant from the 
BGB to Hel. In this case we still use all the above equations 
to calculate the timescales and the luminosity evolution but 
the corresponding value of Mc is a dummy variable. The 
actual core mass at the BGB is given by a mass-dependant 
formula 



Mc.BGB = min 0.95A/c,bagb, (C + ciM" 



(44) 



with C = Mc(Lbgb(Mhop))^ - 
formula is continuous with the 
MhsF, and Mc,bagb given by eq, 

9.20925 10"'^ and C2 = 5.402216 are independent of Z 
that for large enough M we have Mc,bgb ^ 0.098M^'^'^ 
dependent of Z. Thus on the GB we simply take 



ciMhof'^^, ensuring the 
Mc-L relation at M = 
(|6^. The constants ci = 
" " so 



,GB = Mc,BGB + (Mc,HoI - Mc,bgb) T M > MhoF (45) 



with 



t — tsGB 



tv 



to account for the small growth of the non-degenerate core 
while Mc,GB is given by eq. (pST) for M < MueF- A^c.HoI is 



described in Section 5.3 



Furthermore, as giants have a deep convective envelope 
and thus lie close to the Hayashi track, we can find the radius 
as a function of L and M, 



Rot 



where 



A 



(46) 



A = min biM 



and fei ~ 0.4, 62 ~ 0.5 and 63 « 0.7. A useful quantity 
is the exponent x to which 7? depends on M at constant L, 
Rgb (X . Thus we also fit x across the entire mass range 
hy A = bM~^ , ie. a hybrid of bs and 67, to give 

X = 0.30406 + 0.0805C + 0.0897C^ + 0.0878C^ + 0.0222C'' (47) 

so that it can be used if required. Thus for Z = 0.02, as an 
example, we have 

Rgb ~ l.lM-°-^ (L°-^ + 0.383L°'^«) . (48) 

Figure |l^ exhibits the accuracy of eq. ^dj) for solar mass 
models of various metallicity. 



5.3 Core helium burning 

The behaviour of stellar models in the HRD during CHeB 
is fairly complicated and depends strongly on the mass and 
metallicity. For LM stars. He ignites at the top of the GB 
and CHeB corresponds to the horizontal branch (including 
the often observed red clump); the transition between the 
He flash and the start of steady CHeB at the ZAHB is very 
rapid and we take it to be instantaneous. For IM stars, CHeB 
can be roughly divided in two phases, descent along the GB 
to a minimum luminosity, followed by a blue loop excursion 



M = 1-0 
Z - 0.0100 




Figure 12. Relation between radius and luminosity on the giant 
brancli as given by eq. ( ^ ) (solid line) and from the detailed 
models (crosses) for 1.0 Mq at various metallicities. 



to higher Tgff connecting back up to the base of the AGB 
(BAGB). However, not all IM stars exhibit a blue loop, in 
some cases staying close to the GB throughout CHeB (the 
so-called 'failed blue loop'). Sometimes the blue loop is also 
followed by another period of CHeB on the GB but this is 
usually much shorter than the flrst phase and we choose to 
ignore it. For HM stars, He ignites in the HG and CHeB also 
consists of two phases, a blue phase before reaching the GB 
followed by a red (super)giant phase. 

For the purpose of modelling, we define the blue phase 
of CHeB as that part which is not spent on the giant branch. 
This means that the position in the H-R diagram during the 
blue phase can in fact be quite red, e.g. it includes the red 
clump and failed blue loops. By definition, for the LM regime 
the whole of CHeB is blue. For IM stars, the blue phase 
comes after the RG phase, while for HM stars it precedes 
the RG phase. 

The transition between the LM and IM star regime oc- 
curs over a small mass range (a few times 0.1 M0), but it 
can be modelled in a continuous way with a factor of the 
form 1 -I- Qexp 15(M — Mhbf) in the LM formulae (see be- 
low). With a of order unity, this factor can be neglected if 
M <^ MhoF. We also require continuity of LM CHeB stars 
with naked He stars when the envelope mass goes to zero. 
The formulae are also continuous between IM and HM stars 
for Z < 0.002. For higher Z, however, there is a disconti- 
nuity in the CHeB formulae at M = Mfgb, because the 
transition becomes too complicated to model continuously 
while keeping the formulae simple. 

The luminosity at helium ignition is approximated by 



LhcI ~ 



bgM^^" 

1 + ai exp 15(M-Mhof) 

fell +bi2M3-'^ 
613 + M2 



M < A'/hcF 



M > MhcF 



(49) 
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with ai = [fegMHcF^" - LHci(MHcF)]/LHei(MHeF). The ra- 
dius at He ignition is -RhcI = Rgb{M, Lhbi) for M < Mfgb, 
and i?Hoi = ^?mHo for M > max(MFGB, 12.0), with -RmHc 
given by eq. (pa) below. If Mfgb < M < 12.0, we take 



Rn 



f -Rgb(I/Hci) 

\ -RmHc 



log(M/12.0) 

log(MFGB/12.0) 



.(50) 



The minimum luminosity during CHeB for IM stars, 
reached at the start of the blue phase, is given by 



-t/min.Hc — ^Hcl 

with 



bi4 + cM'^^^+°-^ 
6i6 + Mbi5 



(51) 



bi7 



Mfgb" 



+ 



biabi 



Mf 



bis+o.i '' 



so that Z/min,Hc = bn Luci at Af = AfpGB- Continuity with 
HM stars, for which there is no minimum luminosity, is 
achieved by taking bn = 1 for Z < 0.002 (but 617 < 1 for 
Z > 0.002). The radius at this point is -Rgb(M, Lmin.He)- 

For LM stars the ZAHB luminosity Lzahb takes the 
place of Lmin,Hc. To model the ZAHB continuously both 
with the minimum luminosity point at M — MhcF and with 
the naked He star ZAMS (see Section 6.1) for vanishing en- 
velope mass (M — Mc), the ZAHB position must depend on 
Mc as well as M. We define 



M -Mc 



MHeF - Mc ' 

so that < /i < 1, and then take 

izAHB ~ izHc(Mc) + 

1 + &20 bisn' 



(52) 



&19 



02 



l + fe20M^"^^ 1 +a2expl5(M-MHoF)' 

bl8 + LzHajMc) — Lmin,Hc(MHcF) 
Hg(Mhcf) — LzHc{Mc) 



(53) 



where Lzho is defined by eq. (|77[). (Note that this 02 is not 
a constant but depends on M^.) For the ZAHB radius we 
take 



^ZAHB = (1 — /)i?ZHc(Mc) + /i?GB(izAHB); 



(54) 



/ 



(1.0 -I- b2l)M " 

1.0 -f 621^^23 



This formula ensures, apart from continuity at both ends, 
that J?zAHB is always smaller than the GB radius at Lzahb. 

The minimum radius during the blue loop is approxi- 
mated by 



R 



b24M +{b25M)^^0M^^ 



627 -I- M&28 

Then for M < MhoF, we simply take 

i?mHc(AfHcF 



M > A^HoF. 



(55) 



^mHc — ^GB(izAHB) 



-RGB(izAHB(MHcF)) 



to keep iZmHo continuous. 

The luminosity at the base of the AGB (or the end of 
CHeB) is given by 



629M"30 



i/BAGB 



14 

bsi 



a3expl5(M-MHeF) 

+ b32M&33+1.8 



&34 + 



M < Mhof 



M > Mhcf 



(56) 



with aa = [&29MHeF''^° - iBAGB(AW)]/LBAGB(AfHoF). 

The radius at the BAGB is simply i?AGB(Af, Lbagb), as 
given by eq. (|74|). 

The lifetime of CHeB is given by 



bag + {tHcMs(Mc) - 639} (1 
X [1 + 04 exp 15(A//-Mhof)] 
b4 



M)^ 



ifiGB 



biiM"^^ +643M^ 



M < Mu 



M > Mhof 



(57) 



&44 + IVP 

with Q4 = [fHc(AfHoF) — b39]/&39. The term involving 
tHcMs(Mc) ensures continuity with the lifetime of a naked 
He star with M = Mc as the envelope mass vanishes. The 
hfetime of the blue phase of CHeB relative to tne depends in 
a complicated way on M and Z, it is roughly approximated 
by 

M < Mhof 



Tbl ^ 



645 



(1 



M 



Ml 



FGB 



Clbl 



(log- 



M 



/bi(M) 



A/fgb > 
Mhof < M < Mfgb 

M > AfFGB 



(58) 



/bl(MFGB) 

truncated if necessary to give < rbi < 1, where 



and 



645 



/Mhop y 

' \ Mfgb J 



log 



/b,(Af) = M° 



i?mHc(M) 



Mfgb . 

1 big 



-Ragb(Lhci(M)) 



The second term in the IM part of eq. ( p8| ) with a^i as de- 
fined ensures that rbi = 1 at M = MhoF. By taking 645 = 1 
for Z < 0.002 we also have Tbi = 1 at M = Mfgb- The HM 
part also yields Tbi — 1 at M = Mfgb for Z < 0.002, so 
that the transition is continuous for low Z. For Z > 0.002 
the transition is regretably discontinuous. Finally, the ra- 
dius dependence of /bi ensures that rbi = at the same 
mass where iimHc = ^AGB(iHei), i-e. where the blue phase 
vanishes. 

During CHeB, we use the relative age r — (t — tHcO/^Hc 
which takes values between and 1. We define as the 
relative age at the start of the blue phase of CHeB, and Lx 
and Rx are the luminosity and radius at this epoch. Hence, 
Tx = for both the LM and HM regime, and Tx ~ 1 — Tbi 
for IM stars, 




M < MHeF 

Mhsf < M < Mfgb 
M > Mfgb 



and 



Rx — 



RzAHB 

JZGB(imin,Hc) 
RhcI 



M < Mhof 

Mhof < M < Mfgb 

M > Mfgb 



(59) 



(60) 
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Then the luminosity during CHeB is modeUed as 

'isAGB \ ^ ^ 1 

" < r < 1 



(61) 



where 
A = 



1 



e = min(2.5, max(0.4, RmKc/R^)), (62) 



(63) 



The actual minimum radius during CHeB is -Rmin ~ 
mm{RinHc, Rx), because eq. ( ^ ) for -RmHc can give a value 
that is greater than R^ (this property is used, however, to 
compute 4 above). Furthermore, we define Ty as the relative 
age at the end of the blue phase of CHeB, and Ly and Ry as 
the luminosity and radius at r = r^. Hence, Ty = 1 for LM 
and IM stars and Ty = Tbi for IM stars. Ly is given by eq. dHl] ) 
{Ly = I/BAGB for M < Mpgb), and Ry = RAOBiLy). The 
radius during CHeB is modelled as 



R 



where 



Rgb{M,L) 
Ragb{M,L) 

Rmin exp(|p|''') 



0<T<Tx 
Ty<T<l 
Tx < T < Ty 



(64) 



P=i^r.^Y(l^)-(ln^Y(l^). (65) 

V -Kmin / \Ty — TxJ \ Rmin ) \Ty — Tx J 

The core mass Mc,hcI at He ignition is given by the 
Mc-L relation for LM stars, while for M > MhcF the same 
formula can be used as for the BGB core mass (eq. ^) re- 
placing Mc(I/bgb(Mhop)) with A/c(iHoi(MHop)) to ensure 
continuous transition at M = MhcF- For M > 3Mq, McHcI 
is nearly equal to Mc,bgb . The core mass at the BAGB point 
is approximated by 



Mc.bagb = {bseM"^' + 63 



(66) 



where 636 ~ 4.36 x 10"", 637 ~ 5.22 and 638 ~ 6.84 x 10"^. 
In between the core mass is taken to simply increase linearly 
with time 



Mc 



(1 - r)Mc,Hoi + rAfcBAGB. 



(67) 



5.4 Asymptotic giant branch 



During the EAGB, when the H-burning shell is extinct, the 
(H-exhausted) core mass Mc,ho (which we have been calling 
Mc so far because it was the only significant core) stays con- 
stant at the value A/c,bagb. Within the H-exhausted core a 
degenerate carbon-oxygen core, Mc,co, has formed and be- 
gins to grow. At a time corresponding to the second dredge- 
up phase the growing Mc,co catches the H-exhausted core 
and the TPAGB begins. From then on Mc.co and Mc,He are 
equal and grow at the same rate (we neglect the mass, about 
0.01 M0, of the thin helium layer between the two burning 
shells) . 

So on the EAGB we set 

Mc = A'/c,Hc = Mc,BAGB . 

Inside this core, Mc,co grows by He-shell burning, at a rate 
© 1999 RAS, MNRAS 000, 



dictated by the Mc-L relation. Thus we can compute the 
evolution of Mc,co and L in the same way as was done 
for GB stars using eqs. (||) and (||) with Mc replaced by 
Mc,co, tsGB replaced by iBAGB (= tud + tee) and Lbgb 
replaced by Lbagb. We also need to replace An with the 
value appropriate for He burning, Aue- The detailed models 
(Pols et al. 1998) on the EAGB show that the carbon-oxygen 
core is composed of 20% carbon and 80% oxygen by mass so 
for every 4 carbon atoms produced by the triple-a reaction, 
3 will capture an a particle and be converted to oxygen. 
Thus 



E : 



e3a + 0.75ecQ 



15 m{H) 



8.09 X 10"' ere 



so that 



{EXuc 



7.66 X lO"'^ Mq Lo"^Myr"^ 



(68) 



using Xhb ~ 0.98. Although massive stars (M ^ 8) do not 
actually follow a Mc-L relation for the CO core, by making 
the proper (ad hoc) assumptions about the constants in the 
relation, we can still effectively model their evolution in the 
same way as for true AGB stars. 

As already mentioned, the EAGB ends when the the 
growing CO-core reaches the H-exhausted core. If 0.8 < 
A/c,bagb < 2.25, the star will undergo a second dredge- 
up phase at the end of the EAGB phase. During this second 
dredge-up the core mass is reduced to 



Mc,Dv = 0.44Afc,BAGB + 0.448. 



(69) 



We assume that the second dredge-up takes place instanta- 
neously at the moment when Afc,co reaches the value Mc,dv, 
so that also Mc,co = Mc at that point (but note that there 
is then a sudden discontinuity in Mc — Afc.Ho). Similarly, for 
Afc,BAGB < 0.8, the EAGB ends when A/c,co reaches Mc,hc 
without a second dredge-up, ie. Afc,Du ~ Mc^bagb- Stars 
with Afc,BAGB > 2.25 do not undergo second dredge-up, as 
they can ignite carbon non-degenerately, and their evolution 
terminates before they ever reach the TPAGB. 

To determine when the transition from EAGB to 
TPAGB occurs we can simply insert Mc,ov into the Mc-L 
relation to find Ldu. Then we calculate 



tinf,l - (p-IJAhcD (^)^!^ Luv<Lx ^ ^^^^ 



Thus if t > tou the TPAGB has begun and the H-exhausted 
and He-exhausted cores grow together as a common core. 
Once again the Mc-L relation is obeyed and once again we 
can use it in the same way as we did for GB stars if we 
replace tBGB by tr,u and Lbgb by Ldu. As we have both 
hydrogen and helium shell burning in operation then we 
must also replace An by an effective combined rate ^H,He 
where 



A 



AhAi 



H,Hc 



1.27 X 10"^ Mq L0"^Myr"^ . (71) 



+ Anc 

There is however an added complication that it is possible 
for Z/Du > Lx- In this case tinf.i and tx are not needed and 
iinf,2 is given by 

tini,2 = too + 7 — FT ( T— ) ' • (72) 



1 / B \ — 

+ (g-l)AH,HcB \L^) 
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Figure 13. Radius evolution from the ZAMS to the end of the 
AGB for a 5.0 Mq star, for metallicities 0.001 and 0.02, showing 
the detailed model points (solid lines) and the fitted tracks (dash- 
dot lines). 



In this way the L evolution (and thus the R evolution) re- 
mains continuous through the second dredge-up. 

On the TPAGB we do not model the thermal pulses in- 
dividually, but we do take into account the most important 
effect of the thermally pulsing behaviour on the long-term 
evolution, namely that of third dredge-ups. During each in- 
terpulse period, the He core grows steadily, but during the 
thermal pulse itself the convective envelope reaches inwards 
and takes back part of the mass previously eaten up by the 
core. The fraction of this mass is denoted by A. Frost (1997) 
shows that models with 4 < M < 6 and 0.004 < Z < 0.02 
have similar overall behaviour in A where A increases quickly 
and reaches approximately 0.9 after about 5 pulses at which 
it stays nearly constant for the remaining pulses. For lower- 
mass stars there is no evidence for such a high A with a value 
of 0.3 more likely for models of approximately solar mass 
and then a steady increase of A with M to reach Amax ~ 0.9 
before M = 4 (Lattanzio 1989; Karakas et al. 1999) Thus 
we simply take A as constant for each M without any Z 
dependence, 

A = min(0.9,0.3-^0.001M'^). (73) 

Hence, the secular growth of the core mass is reduced with 
respect to that given by the Mc-L relation by a fraction 
A. On the other hand, detailed calculations show that the 
luminosity evolution with time follows the same relation as 
without third dredge- up (Frost 1997), ie. it keeps following 
eqs. ( ^ ) and (^^ as if Mc were not reduced by dredge-up. 
In other words, the Mc-L relation is no longer satisfied in the 
presence of third dredge-up, but we can use it nevertheless to 
compute the evolution of L, while Adc is modified as follows: 

Mc = Mc.DU + (1 - A)(Mc' - Afc.Du), 



where Mc' is from the Mc-L relationship, with no dredge-up, 
and Mc.Du is the value of Mc at the start of the TPAGB. 

The radius evolution is very similar to that of the GB, 
as the stars still have a deep convective envelope, but with 
some slight modifications. The basic formula is the same, 

Ragb ^ A (^L^' + b2L^'-°^ (74) 

where indeed bi and 62 are exactly the same as for 7?gb and 
650 = 65563 for M > MhcF. Also for M > MhcF 

A = min (^651 A/-''=^ , 653M-''='') 

which gives 

Ragb = 1.125M-°-^^ + 0.383L° '"^) , 

as an example, for Z — 0.02. For Al < MhoF the behaviour 
is slightly altered so we take 

650 = 63 
A = 656 + 657M 

for M < A/hcF — 0.2 and linear interpolation between the 
bounding values for MhcF — 0.2 < M < MucF, which means 
that for AI = 1.0 and Z — 0.02 the relation gives 

Ragb ~ 0.95 {L°-'^ + 0.383L°-^^) . 

In Figure |l^ we show the radius evolution of a 5.0 Mq 
star, for Z = 0.001 and Z = 0.02, from the ZAMS to the 
end of the AGB, from both the rapid evolution formulae 
and the detailed models. The AGB phase of the evolution 
is recognised by the sharp increase in radius following the 
phase of decreasing radius during the CHeB blue loop. An 
accurate fit to the AGB radius is required if the formulae 
are to be used in conjunction with binary evolution where 
factors such as Roche-lobe overflow and tidal circularisation 
come into play. In actual fact Fig. |l^ shows that we acheive 
an accurate fit for all phases of the evolution. 

We have now described formulae which cover all phases 
of the evolution covered by the detailed grid of stellar mod- 
els. Figures |l^ and |l^ show synthetic HRDs derived from 
the formulae and are designed to be direct comparisons to 
Figs. |l| and ^ respectively. The excellent performance of the 
fitting formulae is clearly evident. 



6 FINAL STAGES AND REMNANTS 

The AGB evolution is terminated, if not by complete loss of 
the envelope, when the CO-core mass reaches a maximum 
value given by 

Mc,sN = max(A/ch, 0.773Afc,BAGB - 0.35). (75) 

When this maximum core mass is reached before the enve- 
lope is lost, a supernova explosion is assumed to take place. 
For stars with Ai^c.BAGB < 2.25, this should occur during the 
TPAGB phase. In practice mass loss will prevent it from do- 
ing so in most cases of single star evolution, but it may oc- 
cur as a consequence of binary evolution. For such stars, we 
make a further distinction based on whether Afc,BAGB ex- 
ceeds 1.6 M0. For Afc.BAGB < 1.6, when the CO-core mass 
reaches Mch carbon ignites in a degenerate flash, leading 
to a thermonuclear explosion. It is uncertain whether we 
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should expect this to occur for normal SSE but if it does 
then the supernova would be something like "type Ila" (la 
+ hydrogen) and we assume that such a supernova leaves 
no stellar remnant. 

For 1.6 < A/c,BAGB < 2.25, the detailed models show 
that carbon ignites off-centre under semi-degenerate condi- 
tions when Mc,co ^ 1-08 (Pols et al. 1998). Carbon burning 
is expected to lead to the formation of a degenerate ONe- 
core (Nomoto 1984), while the star continues its evolution 
up the AGB. When the core mass reaches Mch, the ONe- 
core collapses owing to electron capture on Mg^'* nuclei. The 
resulting supe rnova explosion leaves a neutron star remnant 
(Section B.2.2). The limiting A^cBAGb values of 1.6 M0 and 
2.25 M0 correspond to initial stellar masses denoted tra- 
ditionally by the symbols Mup and A/oc, respectively. The 
values of A/up and M^c depend on metallicity (see Table 1 
of Pols et al. 1998), this dependence follows from inverting 
eq. ( |66| ) for the values Afc.BAGB = 1-6 and 2.25, respectively. 

If the envelope is lost before Mc reaches Afc.SN (= A/ch) 
on the TPAGB, the remnant core becomes a white dwarf. 
This will be the case for almost all cases of normal SSE. 
For A/c,BAGB < 1.6, this will be a CO white dw arf, fo r 
Afc.BAGB > 1.6 it will be a ONe white dwarf (Section 6.2.1 ). 

Stars with Afc,BAGB > 2.25 develop non-degenerate 
CO-cores which grow only slightly before undergoing central 
carbon burning, rapidly followed by burning of heavier ele- 
ments. Here, Afc sn is the CO-core mass at which this burn- 



ing takes place, because the core mass does not grow signif- 
icantly after C burning. Very quickly, an Fe-core is formed 
which collapses owing to photo-disintegration, resulting in a 
supernova explosion. The supernova leaves either a neutron 
star or, for very massive stars, a black hole (Section |6.2.2 ). 
We assume that a black hole forms if Afc,sN > 7.0, corre- 
sponding to A/c,BAGB > 9.52. 

This means that the lowest mass star to produce a NS 
has an initial mass Af, in the range A^up < M, < Mac with 
the actual value of Af, depending greatly on the mass loss 
rate. Observations would tend to suggest that M, ~ Afcc 
(Elson et al. 1998) and inde ed we find that with our adopted 
mass loss rate (Section 7.1) almost all cases of SSE result in 



WD formation for M < A4c. 

While most stars have their nuclear burning evolution 
terminated on the TPAGB we must make allowances for 
cases of enhanced mass loss, e.g. owing to binary evolution 
processes, that result in termination at an earlier nuclear 
burning stage. If the star loses its envelope during the HG 
or GB phases then the star will become either a HeWD 
(Section 6.2.1), if it has a degenerat e co re {M < A/hcf), 
or a zero-age naked He star (Section |6.l[ ). If during CHeB 
Al = Mc then an evolved naked He star is formed with 
the degree of evolution determined by the amount of central 
helium already burnt. Thus the age of the new star is taken 
to be 
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t ■ 



(76) 



where the primes denote times for the original star and 
t-ReMS is given by eq. (|7^). When the envelope is lost during 
the EAGB so that Aic,Ua = M, a naked helium giant (Sec- 



tion 3.1) is formed as unburnt helium still remains within 
Mc,He through which the growing Mc,co is eating. The age 
of the new star will be fixed by using Mc = M c.co and 
M = Mc,Ho in the HeGB Mc-t relation (see Section O) . We 



note that although naked helium stars are nuclear burning 
stars, ie. not a final state, we still label them as a remnant 
stage because they are the result of mass loss. Also, when 
a WD, NS or BH is formed the age of the star is reset so 
that the remnant b egin s its evolution at zero-age to allow 
for cooling (Section [6. 2| ). 

6.1 Naked helium stars 

The formulae described in this section are based on detailed 
stellar evolution models for naked helium stars, computed 
by one of the authors (ORP) with the same code as used 
for the stellar models described in Section ^. First, a he- 
lium ZAMS of homogeneous models in thermal equilibrium 
was constructed, with composition X — 0, Y = 0.98 and 
Z = 0.02. Starting from this ZAMS, evolution tracks were 
computed for masses between 0.32 and 10 M© spaced by 
approximately 0.1 in log A/. For masses below 2 Mq, the 
tracks were computed until the end of shell He burning and 
for M > 2 Mq , up to or through central carbon burning. 
These models will be discussed in more detail in a forth- 
coming paper (Pols 1999, in preparation). 

The following analytic formulae provide an accurate fit 
to the ZAMS luminosity and radius of naked He stars with 
Z = 0.02: 



LzHa = 



RzHc ~ 



15 262 



Af9 + 29.54Af'^-5 + 31.18M6 + 0.0469 ' 

0.2391Af^-^ 
Af4 -h0.162Af3 +0.0065' 



(77) 



(78) 



The central He-burning lifetime (He MS) is approximated 
by 



^HeMS 



0.4129 + 18.81Ar + 1.853Af'5 



(79) 



The behaviour of L and R during central He burning can be 
approximated by 



I/HeMS = izHe(l -|- 0.45t -|- ttT 

and 

-RhcMS = -RzHe(l + I3t ~ (3t'^). 



(80) 



(81) 



where r — t/tecMS and t is counted from the He ZAMS. a 
and (3 are dependent on mass, as follows: 



a = max(0, 0.85 - 0.08Ar) 
and 

/3 = max(0, 0.4-0.22 log M) . 



(82) 



(83) 



The evolution after the He MS is dominated by the 
growth of the degenerate C-O core for low-mass stars, and 
by evolution up to carbon burning for M ^ 2. Low-mass He 



stars follow an approximate core mass-luminosity relation 
(e.g. Jeffery 1988), and we compute their evolution making 
use of this relation just as we do for GB stars (Section 5.2). 
For massive He stars, although they do not properly follow 
such a relation, an ad hoc Mc-L relation can be used to also 
describe their evolution. The following formula works for the 
whole mass range: 

LhcGb = MW{BM,\DM,^) (84) 

with B = 4.1 X lO"* and D = 5.5 x 10*/(l + 0.4Af''). The first 
term models the 'real' Mc-L relation followed by low-mass 
He stars, while the second, mass-dependent term mimics the 
behaviour for high-mass stars. The evolution of L and Afc 
with time is obtained from eq. (^^ and the equivalents of 
eqs. ( p9[]4^ ) with Ah replaced by Auc as given by eq. (|68[), 
^BGB replaced by tncMS, and Lbgb replaced by Lthc. Lthc 
is the value of L at the end of the He MS, i.e. LhoMS given 
by eq. ( ^o| ) at r = 1. The post-HeMS radius can be approx- 
imated by 



RhcGb = MIN(i?i,i?2), 
/ L \°'^ 

Rl = i?ZHc 7 

V ^THc / 

2 + 



-0.02 



cxp 



/ LtHc 



X = 500- 



Af2.5 ' 



i?2 = 0.08L" 



(85) 
,(86) 

(87) 
(88) 



The first term of Ri models the modest increase in radius 
at low mass and/or L, and the second term the very rapid 
expansion and redward movement in the HRD for M ^ 0.8 
once L is large enough. The star is on what we call the naked 
helium HG (HeHG) if the radius is given by _Ri. The radius 
J?2 mimics the Hayashi track for He stars on the giant branch 
(HeGB). We make the distinction between HeHG and HeGB 
stars only because the latter have deep convective envelopes 
and will therefore respond differently to mass loss. 

The final stages of evolution are equivalent to those of 
normal stars, i.e. as discussed in Section pi but with A/c,bagb 
replaced by the He-star initial mass M in eq. ( [ts] ) as well as 
in the discussion that follows it. If Af < 0.7 M©, the detailed 
models show that shell He burning stops before the whole 
envelope is converted into C and O. We mimic this by letting 
a He star become a CO WD when its core mass reaches the 
value 

ATcmax = min(1.45Ar - 0.31, AT), (89) 
as long as A/c,max < Afc.SN. 

6.2 Stellar remnants 

6.2.1 White dwarfs 

We distinguish between three types of white dwarf, those 
composed of He (formed by complete envelope loss of a 
GB star with M < MhcF, only expected in binaries), 
those composed of C and O (formed by envelope loss of 
a TPAGB star with M < Afup, see above), and those com- 
posed mainly of O and Ne (envelope loss of a TPAGB star 
with Afup < M < Afec). The only distinction we make be- 
tween CO and ONe white dwarfs is in the way they react 
to mass accretion. If A4"wd + Afacc > Mch, after accreting 
an amount of mass Afacc, then a CO WD explodes without 
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leaving a remnant while an ONe WD leaves a neutron star 
remnant with mass A/ns = 1-17 + 0.09 (AfwD + Mace) (see 



later in Section S.2.2). The Chandrasekhar mass is given by 



5-8 , 



so it is composition dependent but the mean molecular 
weight per electron is /ic ~ 2, except for low-mass MS stars 
in cataclysmic variables, so we use Mch = 1-44 at all times. 

The luminosity evolution of white dwarfs is modelled 
using standard cooling theory (Mestel 1952), see Shapiro & 
Teukolsky (1983, pg. 85): 

_ 635MZ°'^ 

^WD — TTT" 1 A . (yU) 



[A(t + 0.1)]'-*' 

where t is the age since formation and A is the effective 
baryon number for the WD composition. For He WDs we 
have = 4 for CO WDs A = 15 and for ONe WDs A = 17. 
Eqn. ( po[ ) is adequate for relatively old WDs. The addition 
of a constant in the factor (t + O.l) mimics the fact that the 
initial cooling is rather faster than given by Mestel theory, 
as well as ensuring that it doesn't start at infinite L, so that 
we effectively start the evolution at a cooling age of 10^ yr. 
Note that the initial cooling of the WD is modelled by the 
small-e nvel ope pertubation functions on the TPAGB (see 
Section 6.3). 



The radius of a white dwarf is given by 



as in Tout et al. (1997). 



6.2.2 Neutron stars and black holes 

When a neutron star or black hole is formed in one of the 
situations given above, we assume that its gravitational mass 
is given by 



Mns = 1.17-|-0.09Mc,sN, 



(92) 



where Mc,sn is the mass of the CO-core at the time of 
supernova explosion. With eq. ([TSj), this leads to a mini- 
mum NS mass of 1.3 Mq , and the criterion for BH formation 
Mc,sN > 7.0 gives a maximum NS mass and minimum BH 
mass of 1.8 Mq. 

The NS cooling curve is approximated by assuming that 
photon emission is the dominant energy loss mechanism, 
which should be true for t > 10^ yrs (see Shapiro & Teukol- 
sky, pg. 330); 

L.s^ 0-02M-^\. (93) 

The upper limit is calibrated to give Toff ^ 2 X HfK which 
is appropriate for the Crab Pulsar and is set constant for 
the first 10^ yrs to refiect the scatter in the observations of 
Tcft for pulsars with an age less than 10^ yrs. Eqn. (^) also 
ensures that Lns < iwD at all times and that neutron stars 
will cool faster than white dwarfs. 

The radius of a NS is simply set to 10 km, i.e. _Rns = 
1.4 10"^ 

We take the black hole radius as the Schwarzschild ra- 
dius 



i?3H = ^^=4.24 10-«Mbh. 



(94) 



The corresponding luminosity of a BH is approximately 
given by 

1.6 X 10"^^° 



(95) 



(Carr & Hawking 1974) which will be negligible except for 
extremely low mass objects and thus we actually set 

-10 



^BH = 10" 



(96) 



to avoid floating point division by zero. 

Note that for all remnants we set Mc = M for conve- 
nience. 



6.3 Small envelope behaviour and hot subdwarfs 

In general the equations in Section ^ accurately describe 
the nuclear burning evolution stages as outlined by our grid 
of detailed models. However, we also find it necessary to 
add some pertubation functions which alter the radius and 
luminosity when the envelope becomes small in mass, in 
order to achieve a smooth transition in the HRD towards 
the position of the remnant. Take, for example, the AGB 
radius where 

7?AGB oc M""= 

so that as M decreases due to mass loss from a stellar wind 
Rags will increase and the star moves further to the red 
in the HRD. In actual fact, as the envelope mass (Menv) 
gets very small, the star becomes bluer and moves across 
the HRD to WD temperatures. In the same way we would 
also expect the luminosity growth rate to decrease until the 
luminosity levels off at some approximately constant value 

for small Menv. 

Thus for any nuclear burning evolution stage where 
there is a well defined core and envelope (i.e. not the MS), 
we define 

where Lo — 7.0 x lO**, k — —0.5 for normal giants and 



Mc, 



Mc 



Mc. 



(98) 



for helium giants. Then if < 1.0 we perturb the luminosity 
and radius using 

L 



L' = Lc 



Lc 



where 



(1 + b^) (M/fc)^ 
1 + iti/hf 

1 -f (m/c)3 



with 



b = 0.002 max 



(99) 
(100) 

(101) 
(102) 

(103) 
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c = 0.006 max 



q = 



Mi) 



V M 



(104) 
(105) 



The luminosity and radius of the star are then given by L' 
and R'. 

In the above formulae, Lc and Rc are the luminosity 
and radius of the remnant that the star would become if it 
lost all of its envelope immediately. Thus we set M = Mc in 
the appropriate remnant formulae. If the star is on the HG 
or GB then we have, for AI < MhoF, 

Lc = LzHc (Mc) 
Rc = RzHc (Mc) 
otherwise 

Lc = I/WD (Mc) , i.e. eq. (^ with A = 4 and t = 0, 
Rc = Rwo (Mc) . 

During CHeB the remnant will be an evolved helium MS star 
so we use Mc and r = (t — tHci)/tHc in eqns. ( po|) and (jsi]) 
to give Lc and -Rc respectively. On the EAGB the remnant 
will be a helium HG or GB star with M = Mc^hc so that Lc 
comes from the HeGB Mc-L relation with Mc — Mc,co and 
Rc from i?HcGB = (i\fc,Ho,Lc). For the TPAGB, HeHG and 
HeGB the remnant will most likely be a CO WD so 



eq. (^ with A = 15 and f = 0, 



Lc = I/WD (Mc) , 

Rc = Rwo (Mc) . 

Figure ^ shows how a model incorporating mass loss (us- 
ing the prescription outlined in Section 7.1) and the small- 
envelope pertubation functions deviates from a model with- 
out either. No difference is evident until the stellar wind 
becomes appreciable as the star evolves up the AGB. As the 
envelope mass is reduced the star initially moves to the right 
of the AGB becoming redder in accordance with eq. (|74|). 
Then as the envelope is reduced even further in mass the 
star moves to the left in the HRD, under the influence of 
the pertubation functions, becoming bluer as the hot core 
starts to become visible. Thus we have in effect mimicked 
the planetary nebulae nucleus phase of evolution which fin- 
ishes when the star joins up with the white dwarf cooling 
track (marked by a cross on the figure) . The behaviour of the 
core-mass-luminosity relation for the same models is shown 
in Fig. Both the helium and the carbon-oxygen cores are 
shown on the AGB until second dredge-up when the helium 
core is reduced in mass and the two grow together. It can 
be seen that after second dredge-up the slope of the relation 
changes as a result of third dredge-up during the TPAGB 
phase. 

We should note that Rc can be used directly as a fairly 
accurate estimate of the current core radius of the star ex- 
cept when Rc is given by Rwo- In that case nuclear burning 
will be taking place in a thin shell separating the giant core 
from the envelope so that the core will be a hot subdwarf 
for which we assume the radius Rc — 5Rwd (Mc). It is also 
necessary to check that Rc < R in all cases. 



7 MASS LOSS AND ROTATION 
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Figure 16. Synthetic evolution tracks on the HRD for a 5.0 Mq 
star without mass loss (dash-dot line) and with mass loss (points). 
The cross marks where the WD cooling track begins. 
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Figure 17. Relation between core mass and luminosity for a 
5.0 Mq star as given by the formulae without mass loss (dash-dot 
line) and with mass loss (points). Both the helium and carbon- 
oxygen cores are shown for the EAGB phase. 



7.1 Mass loss 

We now describe a particular mass loss prescription which is 
independent of the previous formulae and fits observations 
well. On the GB and beyond, we apply mass loss to the 
envelope according to the formula of Kudritzki & Reimers 
(1978), 
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r p 

Mr = 774 X 10-^3— MQyr~\ (106) 

with a value ol 77 = 0.5. Our value for rj is within the lim- 
its set by observations of Horizontal Branch morphology in 
Galactic globular clusters (Iben & Re nzin i 1983) and we 
don't include a .Z-dependence in eq. (106) as there is no 
strong evidence that it is necessary (Iben & Renzini 1983; 
Carraro et al. 1996). On the AGB, we apply the formulation 
of Vassihadis & Wood (1993), 



log Mvw 



-11.4 + 0.0125[Po - 100 max(M - 2.5, 0.0)], 



to give the observed rapid exponential increase in M with 
period before the onset of the the superwind phase. The 
steady superwind phase is then modelled by applying a max- 
imum of Mvw = 1.36 X 10"''L Moyr"^ Po is the Mira 
pulsation period given by 

log Po = min (3.3, -2.07 - 0.9 log M -^1.94 log R) . 

For massive stars we model mass loss over the entire 
HRD using the prescription given by Nieuwenhuijzen & de 
Jager (1990), 

1/2 



Mn,7 = 



9.6 X 10"""^ 



15 r>0-81 r 1-24 j^0.16 



for L > 4000 L©, modified by the factor Z'^^^ (Kudritzki et 
al. 1989). 

For small H-envelope mass, fi < 1.0, we also include 
a Wolf-Rayet-like mass loss (Hamann, Koesterke & Wes- 
solowski 1995; Hamann & Koesterke 1998) which we have 
reduced to give 

MwR = 10"'^L'-^(1.0-m) Moyr"' 

where /i is given by eqn. (p7|). The reduction is necessary in 
order to produce sufficient black holes to match the number 
observed in binaries. 

We than take the mass loss rate as the dominant mech- 
anism at that time 

M = max (Mb., Mvw, Mnj, Mwr) M0yr"\ 

In addition we add a LBV-like mass loss for stars be- 
yond the Humphreys-Davidson limit (Humphreys & David- 
son 1994), 

Mlbv = 0.1 {W-"RL'^^Y {-^^T^ - ^-o) ^eyr^'' 

if L > 6x10^ and 10"^i?L^/^ > 1.0, so that M = M+Mlbv- 
For naked helium stars we include the Wolf-Rayet-like 
mass loss rate to give 



M = max(MR,Mv 



(m = o)) 



M 



©yr 



The introduction of mass loss means that we now have 
two mass variables, the initial mass Mo and the current mass 
Mt (= M). From tests with mass loss on detailed evolution 
models we found that the luminosity and timescales remain 
virtually unchanged when mass loss is included, during the 
GB and beyond, but that the radius behaviour is very sensi- 
tive. Thus we use Afo in all formulae that involve the calcu- 
lation of timescales, luminosity or core mass and we use Mt 
in all radius formulae. When a MS star loses mass, which 
may occur in a stellar wind for massive stars or as a result 
of mass transfer, it will evolve down along the MS to lower 



L and T^b because of the decrease in central density and 
temperature. The luminosity responds to changes in mass 
because the size of the core depends on the mass of the star 
and therefore Mo, which is more correctly the effective ini- 
tial mass, is kept equal to the current mass while the star 
is on the MS. We must effectively age the star, so that the 
fraction of MS lifetime remains unchanged, by using 

this 

where primes denote quantities after a small amount of mass 
loss {tMs' > tus thus t' > t). Even though the star has been 
aged relative to stars of its new mass, its remaining MS life- 
time has been increased. Naked helium main-sequence stars 
must also be treated in the same way with tus replaced 
by tuoMS. During the giant phases of evolution the age de- 
termines the core mass which will be unaffected by mass 
changes at the surface, as the core and envelope are effec- 
tively decoupled in terms of the stellar structure, so that 
the age and the initial mass do not need to be altered. HG 
stars will respond to changes in mass on a thermal timescale 
and thus, as our detailed models show is necessary, we keep 
Mo — Mt during the HG and the star is aged according to 



t' 



^Ms') 



^MS) 



{t — tyis) 



whenever mass is lost. However, as the core mass depends 
on Mq, see eqs. ( p^ - |3C| ), there exists a limiting value beyond 
which Mq cannot be decreased. To do otherwise would lead 
to an unphysical decrease in the core mass. Therefore our 
treatment of mass loss on the HG is a mixture of the way 
the MS and giant phases are treated which in a sense reflects 
the transitional nature of the HG phase of evolution. 

When a LM star experiences the He-flash and moves to 
the ZAHB we reset Mq = Mt, so that t = tnci (Mq) as it 
is now a new star with no knowledge of its history. We also 
reset Mo = Mt when naked helium star evolution is begun. 



7.1.1 The white dwarf initial-final mass relation 

If a star is to evolve to become a WD the minimum mass pos- 
sible for the WD is the core mass at the start of the TPAGB. 
Thus an accurate empirical relation between white dwarf 
masses and the initial mass of their progenitors provides an 
important calibration of the mas s los s required on the AGB. 
This helps to constrain rj in eq. (106) which, for now, is ba- 
sically a free parameter. The commonly used method to ob- 
tain the initial-flnal mass relation (IFMR) for white dwarfs 
is to use WDs that are members of clusters with known ages. 
Their radii, masses and cooling times can be obtained spec- 
troscopically so that by subtracting the cooling time from 
the cluster age the time spent by the progenitor from the 
ZAMS to the AGB can be estimated. The initial progenitor 
mass. Mi, must then be derived using appropriate stellar 
models so that this a semi-empirical method for deflning 
the IFMR. Using data from WDs in galactic open clusters 
Weidemann (1987) derived such a semi-empirical IFMR as 
shown in Fig. As Jeffries (1997) rightly points out, an 
IFMR derived by this method will be sensitive to the amount 
of core overshooting included in the stellar evolution models. 
The effect of increased overshooting is to decrease the de- 
rived cluster age, thus increasing the progenitor lifetime and 
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Figure 18. Relation between white dwarf mass and the ZAMS 
mass of its' progenitor, i.e. the initial— final mass relation (IFMR). 
The IFMR from the evolution formulae (solid line) is given for 
Z = 0.02 and Z = 0.004 as well as the corresponding core masses 
at the start of the TPAGB (dotted lines). The vertical lines corre- 
spond to Mup. Weidemann's (1987) semi-empirical IFMR (dash- 
dot line) and the NGC 2516 white dwarfs (crosses) from Jeffries 
(1997) are shown. 



decreasing Mi. The IFMR will also be sensitive to changes 
in metallicity. 

Jeffries (1997) presents initial and final masses for 4 
WDs found in the young open cluster NGC 2516 which has 
a metallicity of Z ~ 0.009. The initial progenitor masses 
are derived from the stellar models of Schaerer et al. (1993) 
with Z = 0.008 and moderate core overshooting. We show 
the data points for these 4 WDs in Fig. |l^ as well as the 
IFMR given by our formulae for Z = 0.02 and Z = 0.004 
(the IFMR for Z = 0.009 wiU he between these two), and 
the corresponding core mass at the start of the TPAGB. As 
the TPAGB core mass is the minimum possible mass for the 
WD it is clear that our formulae are in disagreement with the 
semi-empirical IFMR of Weidemann (1987). Jeffries (1997) 
was in similar disagreement with the semi-empirical IFMR. 
However the IFMRs from our formulae are in good agree- 
ment with the NGC 2516 data, taking the associated errors 
of the data points into account. Thus there is no contradic- 
tion with the mass loss prescription used for the formulae 
however, we note that an empirical IFMR is required before 
concrete conclusions can be drawn. 



T.2 Rotation 

As we plan to use the evolution routines for single stars in 
binary star applications it is desirable to follow the evolu- 
tion of the stars' angular momentum. To do this we must 
start each star with some realistic spin on the ZAMS. A 
reasonable fit to the Urot MS data of Lang (1992) is given by 



330M"" , _i 
"-^(^^^ 15.0 + M3-45 km^ (107) 

so that 

n = 45.35 yr-i . (108) 

-RzAMS 

The angular momentum is then given by 
Jspin ^m = kMR^n 

where the constant k depends on the internal structure, e.g. 
fc = 2/5 for a solid sphere and k — 2/3 for a spherical shell. 
In actual fact we find the angular momentum by splitting the 
star into two parts, consisting of the core and the envelope, 
so that 

Jspin = (fe (M - Mc) + fcgA/cEc^) !^ (109) 

where ^2 = 0.1, based on detailed giant models which reveal 
k — O.lMonv/A/, and ks = 0.21 for an n = 3/2 polytrope 
such as a WD, NS or dense convective core. This works well 
for post-MS stars which have developed a dense core whose 
rotation is likely to have decoupled from the envelope while 
also representing the near uniform rotation of homogenous 
MS stars which have Mc = 0.0. When the star loses mass 
in a stellar wind the wind will carry off angular momentum 
from the star at a rate given by 

j = kMh 

where h = R^Q. Thus 

Jspin = Jspin - '^AMR^Q (110) 

when the star loses an amount of mass AAf, where we take 
fc = 2/3 as we assume that all the mass is lost uniformly at 
the surface of the star, ie. from a spherical shell. 

We also include magnetic braking for stars that have 
appreciable convective envelopes where 

jn.b = 5.83 X 10-''^ {Rnf M0R|yr-^ (111) 

with Q in units of years. However, following Rappaport et al. 
(1983), we don't allow magnetic braking for fully convective 
stars, M < 0.35. 

For most stars Menv is simply given by M — Mc however 
the case is slightly more complicated for MS and HG stars. 
Our detailed models show that MS stars are fully convec- 
tive for M < 0.35 so that Mcnv.o = M and that MS stars 
with M > 1.25 have little or no convective envelope so that 
AJonv.o = 0.0, independent of Z. In between we take 

/I 25 — Af \ ^ 
Menv,o = 0.35 i—^ J 0.35 < AJ < 1.25 . 

The convective envelope, if it is present, will diminish as the 
star evolves across the MS so we take 

AJenv = Afenv.O (1.0 - t)^^^ , 

where 
_ t 

and Mcnv,o is effectively the ZAMS value. On the HG we 
assume that the convective core gradually establishes itself 
so that 

Mcnv =t(M - Mc) 
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where 



t — tms 



8 DISCUSSION 

The possible paths of evolution through the various phases 
described in the preceeding sections are illustrated in Fig. ^[ 
In Fig. ^ we show the distribution of remnant masses and 
types, as a function of initial stellar mass, for Population I 
and II stars as given by the rapid evolution code. The dis- 
tribution approximates what we would see if a population 
of single stars were to be evolved to the current age of the 
Galaxy. The variation in behaviour produced by a change in 
metallicity should once again be noted. These variations are 
due to changes in the evolution rates as a function of initial 
mass, brought about by changes in the composition. The ini- 
tial mass above which stars will become black holes rather 
than neutron stars is not well constrained which is why we 
use the maximum AGB core mass in the formulae to decide 
the outcome, corresponding to a transition at Mq ~ 30 Mm 
(varying with metallicity). It can also be seen from Fig. M 
that, above this mass, a small pocket of neutron star for- 
mation occurs in what would normally be assumed to be 
a region of black hole formation on the diagram. This be- 
haviour corresponds to a massive star losing its envelope on 
the HG so that the star enters the naked helium MS phase, 
where the mass loss rate increases, causing a reduction in 
Mq. As a result a lower value than otherwise expected for 
Mns is given by eq. ( p^ ) when the naked helium evolution 
ends. 

The formulae described in this paper are available in 
convenient subroutine form as a SSE package, which we also 
term 'the rapid evolution code', that contains; 

EVOLVE The main routine which, amongst other things, 
initialises the star, chooses the timesteps and implements 
mass loss. 

ZCNSTS Subroutine which sets all the constants of the 
formulae which depend on metallicity so that there is no Z 
dependence elsewhere. This needs to be called each time Z 
is changed. 

STAR Subroutine which derives the landmark timescales 
and luminosities that divide the various evolution stages. 
It also calculates tN which is an estimate of the end of the 
nuclear evolution, ie. when Mc = min (Mt, Mc,sn), assuming 
no further mass loss. 

HRDIAG Subroutine to decide which evolution stage the 
star is currently at and then to calculate the appropriate L, 
R and Mc. 

ZFUNCS Contains all the detailed evolution formulae as 
functions. 

MLWIND derives the mass loss as a function of evolution 
stage and the current stellar properties. 

In the absence of mass loss STAR is only required at the be- 
ginning of the evolution and then HRDIAG can be called at 
any age to return the correct stellar quantities. When mass 
loss is included, HRDIAG must be called often enough that 
only a small amount of mass is lost during each timestep. 
STAR also needs to be called often as some timescales need 
to be reset after changes of type, e.g. start of the HeMS, as 



do some luminosities, e.g. Lzahb depends on the envelope 
mass at the He-flash. 

The following timesteps, 5tk, are assigned according to 
the stellar type, k: 



Stk 



jOO^MS 


fc = 0, 1 


20 (^BGB — tMs) 


fc = 2 


^ (tinf.l - t) 


k — 3 t tx 


5^ (tinf,2 - t) 


k = 3t>U 




fc = 4 


< 50 (^inf.l - t) 


k — 5,6t < tx 


(tinf,2 - t) 


k = 5,6t>U 


2oteGMS 


fc = 7 


^ (tinf.l - t) 


k^8,9t<tx 




k ^ 8,9 t>tx 


max(0.1, lO.Ot) 


fc > 10 



In addition we impose a maximum TPAGB timestep of 
5x 10"'^ Myr so that important contributions from the small- 
envelope pertubation functions are not missed. We also cal- 
culate Ste, the time to the next change of stellar type (e.g. 
5te — t-MS — t for k = 0,1), and St-^ which is the current 
remaining nuclear lifetime of the star (i.e. Stf^ = t-^ — t as- 
suming that the star is in a nuclear burning stage, otherwise 
tN is set to some large dummy value). If necessary we limit 
the timestep such that mass loss will be less than 1% over 
the timestep. 



5tn 



-0.01 



M 



and we also limit the timestep so that the radius will not 
change by more than 10%, 

StR = 0.1-5- . 

\R\ 

Therefore the timestep is given by 

St — min{Stk,Ste,St!^,5tmi,5tR) . (112) 

In some cases the choice of timesteps is purely for aesthetic 
purposes so the size could easily be increased with no loss 
of accuracy if extra speed is required, such as for evolving 
large stellar populations. For example, the MS can be safely 
done in one timestep but then, for an individual star, the 
hook feature would not appear on a HRD plotted from the 
resulting output. 

Using the SSE package we can evolve 10000 stars up to 
the age of the Galaxy in approximately 100 s of cpu time 
on a Sun SparcUltralO workstation (containing a 300 MHz 
processor). Thus a million stars can be evolved in roughly 
the time taken to compute one detailed model track. This 
speed coupled with the accuracy of the formulae make the 
SSE package ideal for any project that requires information 
derived from the evolution of a large number of stars. How- 
ever, the formulae do not render the model grid of Pols et 
al. (1998) redundant as it contains a wealth of information 
detailing the interior structure of each star, information that 
the formulae simply cannot provide. In actual fact the two 
approaches complement one another. 

The evolution formulae described in this paper have 
been incorporated into a rapid binary evolution algorithm 
so that we can conduct population synthesis involving sin- 
gle stars and binaries. The SSE subroutines have also been 
added to an A^-body code for the simulation of cluster pop- 
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Figure 19. 
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Figure 20. Distribution of remnant masses and types after 1.2 X 
10^" yrs of evolution, as a function of initial mass, for Z = 0.001 
(hollow symbols) and Z = 0.02 (filled symbols). 

ulations. In the future we plan to make 5ov a free param- 
eter as a variable amount of convective overshooting may 
be preferable, especially in the mass range of 1.0 to 2.0 Mq. 
Formulae that describe surface element abundances will also 
be added so that the rapid evolution code can be used for 
nucleosynthesis calculations. 

To obtain a copy of the SSE package described in this 
paper send a request to the authors who will provide the 
FORTRAN subroutines by ftp. 
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APPENDIX 

The Z dependence of the coefficients a„ and b„ is given here. Unless otherwise stated 

a„ = a + /3C + 7C^ + vC' + K* , 
and similarly for b^, where 

C = log(Z/0.02) . 
The variables 
a = log(Z) 
and 

P = C + i-0 

are also used. 



ai 1.593890(+3) 

a2 2.706708(+3) 

as 1.466143(+2) 

at 4.141960(-2) 

as 3.426349(-l) 



2.053038(+3) 
1.483131(+3) 
-1.048442(4-2) 
4.564888(-2) 



1.231226(+3) 
5.772723(+2) 
-6.795374(+l) 
2.958542(-2) 



2.327785(+2) 
7.411230(+1) 
-1.391127(+1) 
5.571483(-3) 



ae 1.949814(+1) 

a-r 4.903830(4-0) 

as 5.212154(-2) 

ag 1.312179(+0) 

aio 8.073972(-l) 



1.758178(4-0) -6.008212(4-0) -4.470533(+0) 



3.166411(-2) 
-3.294936(-l) 



-2.750074(-3) 
9.231860(-2) 



-2.271549(-3) 
2.610989(-2) 



a'li 1.031538(+0) 

a'i2 1.043715(4-0) 

ai3 7.859573(4-2) 

ai4 3.858911(4-3) 

ais 2.888720(+2) 

ai6 7.196580(+0) 



-2.434480(-l) 
-1.577474(-H0) 
-8.542048(-H0) 
2.459681(4-3) 
2.952979(-f2) 
5.613746(-1) 



7.732821(+0) 
-5.168234(+0) 
-2.642511(+1) 
-7.630093(4-1) 
1.850341(4-2) 
3.805871(-1) 



6.460705(+0) 
-5.596506(-H0) 
-9.585707(4-0) 
-3.486057(-h2) 
3.797254(4-1) 
8.398728(-2) 



1.374484(-H0) 
-1.299394(+0) 

-4.861703(4-1) 



(2ii — a3^]^ai4 

0-12 — 0,12^14 





a 


P 


7 


■1 






2.187715(-1) 


-2.154437(4-0) 


-3.768678(4-0) 


-1.975518(4-0) 


-3.021475(-1) 




1.466440(4-0) 


1.839725(4-0) 


6.442199(4-0) 


4.023635(4-0) 


6.957529(-l) 


0,20 


2.652091(4-1) 


8.178458(-fl) 


1.156058(4-2) 


7.633811(+1) 


1.950698(+1) 


a,2i 


1.472103(4-0) 


-2.947609(-H0) 


-3.312828(4-0) 


-9.945065(-l) 




(122 


3.071048(4-0) 


-5.679941 (-H0) 


-9.745523(4-0) 


-3.594543(4-0) 




023 


2.617890(4-0) 


1.019135(4-0) 


-3.292551(-2) 


-7.445123(-2) 




024 


1.075567(-2) 


1.773287(-2) 


9.610479(-3) 


1.732469(-3) 




125 


1.476246(4-0) 


1.899331(+0) 


1.195010(-F0) 


3.035051(-1) 




126 


5.502535(4-0) 


-6.601663(-2) 


9.968707(-2) 


3.599801(-2) 





logaiv = max (0.097 - 0.1072 ((7 4- 3), max (0.097, min (0.1461, 0.1461 4- 0.1237 (a 2)))) 

0-19 = Cil9'^20 
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a29 



ass 
ass 





a 


P 


7 


n 




9.511033(+1) 


6.819618(+1) 


-1.045625(+1) 


-1.474939(+1) 


^28 


3.113458(+1) 


1.012033(+1) 


-4.650511(+0) 


-2.463185(+0) 


/ 

0^29 


1.413057(+0) 


4.578814(-1) 


-6.850581(-2) 


-5.588658(-2) 


aso 


3.910862(+1) 


5.196646(+1) 


2.264970(+l) 


2.873680(+0) 


asi 


4.597479(+0) 


-2.855179(-1) 


2.709724(-l) 




as2 


6.682518(+0) 


2.827718(-1) 


-7.294429(-2) 




_ 132 
— "29 






P 


7 




as4 


1.910302(-1) 


1.158624(-1) 


3.348990(-2) 


2.599706(-3) 




3.931056(-1) 


7.277637(-2) 


-1.366593(-1) 


-4.508946(-2) 


asG 


3.267776(-l) 


1. 204424 (-1) 


9.988332(-2) 


2.455361(-2) 


asi 


5.990212(-1) 


5.570264(-2) 


6.207626(-2) 


1.777283(-2) 


— inin 


(1.4,1.5135 + 0.37690 






= max 


; (0.6355 - 0.4192C, max (1.25, 033)) 






a 


P 


7 




ass 


7.330122(-1) 


5.192827(-1) 


2.316416(-1) 


8.346941(-3) 


as9 


1.172768(+0) 


-1.209262(-1) 


-1.193023(-1) 


-2.859837(-2) 




3.982622(-l) 


-2.296279(-l) 


-2.262539(-l) 


-5.219837(-2) 


an 


3.571038(+0) 


-2.223625(-2) 


-2.611794(-2) 


-6.359648(-3) 


a42 


1.9848(+0) 


1.1386(+0) 


3.5640(-l) 




O43 


6.300(-2) 


4.810(-2) 


9.840(-3) 




a44 


1.200(+0) 


2.450(+0) 






= mill 


(1.25, max (1.1, 


CI42)) 






= min 


(1.3, max (0.45, 


CI44)) 








a 


P 


7 




045 


2.321400(-1) 


1.828075(-3) 


-2.232007(-2) 


-3.378734(-3) 




1.163659(-2) 


3.427682(-3) 


1.421393(-3) 


-3.710666(-3) 


a47 


1.048020(-2) 


-1.231921(-2) 


-1.686860(-2) 


-4.234354(-3) 


a4s 


1.555590(+0) 


-3.223927(-l) 


-5.197429(-1) 


-1.066441(-1) 


^49 


9.7700(-2) 


-2.3100(-1) 


-7.5300(-2) 




a^o 


2.4000(-l) 


1.8000(-1) 


5.9500(-l) 




asi 


3.3000(-l) 


1.3200(-1) 


2.1800(-1) 




O52 


1.1064(+0) 


4.1500(-1) 


1.8000(-1) 




053 


1.1900(+0) 


3.7700(-l) 


1.7600(-1) 





042 

^44 



a49 = max (a49, 0.145) 

flso = min (aso, 0.306 + 0.053C) 

flsi = min (asi, 0.3625 + 0.062C) 

052 = max (a52, 0.9) 

^52 = min (0152, 1.0) for Z > 0.01 

ass = max [a^s, 1-0) 

as3 = min (as3, 1.1) for Z > 0.01 
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a 


P 


7 


V 






3.855707(-l) 


-6.104166(-1) 


5.676742(4-0) 


1.060894(+1) 


5.284014(+0) 




3.579064(-l) 


-6.442936(-l) 


5.494644(+0) 


1.054952(+1) 


5.280991 (+0) 




9.587587(-l) 


8. 777464 (-1) 


2.017321(-1) 






157 


1.5135(+0) 


3.7690(-l) 









= mill (1.4, 057) 
057 = max (0.6355 - 0.4192C, max (1.25, 057)) 





a 


P 


7 




ass 


4.907546(-l) 


-1.683928(-1) 


-3.108742(-1) 


-7.202918(-2) 




4.537070(+0) 


-4.465455(+0) 


-1.612690(+0) 


-1.623246(+0) 


1^60 


1.796220(+0) 


2.814020(-1) 


1.423325(+0) 


3.421036(-1) 


161 


2.256216(+0) 


3. 773400 (-1) 


1.537867(+0) 


4.396373(-l) 


062 


8.4300(-2) 


-4.7500(-2) 


-3.5200(-2) 




^63 


7.3600(-2) 


7.4900(-2) 


4.4260(-2) 






1.3600(-1) 


3.5200(-2) 






^65 


1.564231(-3) 


1.653042(-3) 


-4.439786(-3) 


-4.951011(-3) -1.216530d-03 




1.4770(+0) 


2.9600(-l) 






^67 


5.210157(+0) 


-4.143695(+0) 


-2.120870(4-0) 




^68 


1.1160(+0) 


1.6600(-1) 







= max (0.065, 002) 

063 = min (0.055, O63) for Z < 0.004 

064 ~ max (0.091, min (0.121, 064)) 

066 ~ max (066, min (1.6, —0.308 — 1.046^)) 
O66 ~ max (0.8, min (0.8 — 2.0^, O66)) 
068 = max (0.9, min (oes, 1.0)) 
a64 — B — a a (M — O66) for O68 > Oee 
O68 — min (0,38,066) 



P 



069 
O70 
O71 
O72 
O73 
O74 



1.071489(+0) 
7.108492(-1) 
3.478514(+0) 
9.132108(-1) 
3.969331(-3) 
1.600(+0) 



-1.164852(-1) 
7.935927(-l) 
-2.585474(-2) 
-1.653695(-1) 
4.539076(-3) 
7.640(-l) 



-8.623831(-2) 
3.926983(-l) 
-1.512955(-2) 

1.720906(-3) 
3.322(-l) 



-1.582349(-2) 
3.622146(-2) 
-2.833691(-3) 
3.636784(-2) 
1.897857(-4) 



O72 = max (072, 0.95) for Z > 0.01 
074 = max (1.4, min (074, 1.6)) 





a 


P 


7 


V 


075 


8.109(-1) 


-6.282(-l) 






076 


1.192334(-2) 


1.083057(-2) 


1.230969(-K0) 


1.551656(-f0) 


077 


-1.668868(-1) 


5.818123(-1) 


-1.105027(+1) 


-1.668070(4-1) 


078 


7.615495(-1) 


1.068243(-1) 


-2.011333(-1) 


-9.371415(-2) 


O79 


9.409838(4-0) 


1.522928(4-0) 






O80 


-2.7110(-1) 


-5.7560(-l) 


-8.3800(-2) 




081 


2.4930(4-0) 


1.1475(4-0) 
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flys ~ max (1.0, min (075, 1.27)) 

075 = max (075, 0.6355 - 0.4192^) 

076 = max [a-re, -0.1015564 - 0.2161264C - 0.05182516C^) 

077 = max (-0.3868776 - 0.5457078C - 0.1463472C^ min (0.0, 077)) 

078 = max (0.0, min (075, 7.454 + 9.046C)) 

079 = min (a79, max (2.0, —13.3 — 18. 6(^)) 
Ago ~ max (0.0585542, asa) 

flgi — min (1.5, max (0.4, Osi)) 





a 


/3 


7 


■n 






3.9700(-l) 


2.8826(-l) 


5.2930(-l) 






hi 


9.960283(-l) 


8.164393(-1) 


2.383830(+0) 


2.223436(+0) 


8.638115(-1) 


bs 


2.561062(-1) 


7.072646(-2) 


-5.444596(-2) 


-5.798167(-2) 


-1.349129(-2) 


be 


1.157338(+0) 


1.467883(+0) 


4.299661(+0) 


3.130500(+0) 


6.992080(-l) 


br 


4.022765(-l) 


3.050010(-1) 


9.962137(-1) 


7.914079(-1) 


1.728098(-1) 



61 = min (0.54, 61) 

_ -1^0-4. 6739-0.93940- 

62 = min (max (62, -0.04167 + 55.67Z) , 0.4771 - 9329.21^^-^'') 

63 = max (-0.1451, -2.2794 - 1.5175cr - 0.254(7^) 
63 = 10^3 

63 = max (63,0.7307 + 14265. l^^'^^^) for Z > 0.004 

64 = 64 -f 0.1231572C^ 
66 = 66 + 0.01640687C'^ 





Q 




7 


69 


2.751631(+3) 


3.557098(+2) 




610 


-3.820831(-2) 


5.872664(-2) 






1.071738(4-2) 


-8.970339(+l) 


-3.949739(+l) 


612 


7.348793(+2) 


-1.531020(+2) 


-3.793700(+l) 


6'l3 


9.219293(+0) 


-2.005865(+0) 


-5.561309(-1) 



bii = 6'ii^ 
613 = 6'i3^ 



Q /3 7 

614 2.917412(+0) 1.575290(+0) 5.751814(-1) 

615 3.629118(+0) -9.112722(-1) 1.042291(+0) 
6'i6 4.916389(+0) 2.862149(+0) 7.844850(-l) 



614 = b'J'' 

616 = b'J'^ 

617 = 1.0 

617 = 1.0 - 0.3880523 (C + l.O)^'®'^^^'*^ for C > -1.0 
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a 


P 


7 


V 


bi8 


5.496045(+l) 


-1.289968(4-1) 


6.385758(+0) 






1.832694(+0) 


-5.766608(-2) 


5.696128(-2) 




620 


1.211104(+2) 








02 1 


01 /inoo/ 1 o^ 
Z.214Uoo(+2) 


10^7110/ 1 o\ 
2.18/ llo(42j 


1.1 (Ul ( 1 (+1 ) 


-2.03o340(41j 


622 


2.063983(+0) 


7.363827(-l) 


2.654323(-l) 


-6.140719(-2) 


623 


2.003160(+0) 


9.388871(-1) 


9.656450(-l) 


2.362266(-l) 


b'24 


1.609901(4-1) 


7.391573(40) 


2.277010(+1) 


8.334227(+0) 


625 


1. 747500 (-1) 


6.271202(-2) 


-2.324229(-2) 


-1.844559(-2) 


627 


2.752869(+0) 


2.729201(-2) 


4.996927(-l) 


2.496551(-1) 


628 


3.518506(+0) 


1.112440(40) 


-4.556216(-1) 


-2.179426(-1) 



b24 = b'J'^ 

feae = 5.0 - 0.09138012Z-' 
627 = 6^7'^'' 





Q 


P 


7 


V 


&29 


1.626062(42) 


-1.168838(+1) 


-5.498343(+0) 




630 


3.336833(-l) 


-1.458043(-1) 


-2.011751(-2) 




631 


7.425137(41) 


1.790236(41) 


3.033910(+1) 


1.018259(+1) 


632 


9.268325(42) 


-9.739859(+l) 


-7.702152(+1) 


-3.158268(41) 


633 


2.474401(40) 


3.892972(-l) 






^'34 


1.127018(41) 


1.622158(40) 


-1.443664(+0) 


-9.474699(-l) 



b3i = b'J^' 
634 = b'J'^ 



a f3 J r/ 

636 1.445216(-1) -6.180219(-2) 3.093878(-2) 1.567090(-2) 

637 1.304129(40) 1.395919(-1) 4.142455(-3) -9.732503(-3) 

638 5.114149(-1) -1.160850(-2) 



^36 = &36 
637 = 4.O637 

bas = fess 





a 


P 


7 




&39 


1.314955(42) 


2.009258(41) 


-5.143082(-1) 


-1.379140(40) 


640 


1.823973(41) 


-3.074559(-H0) 


-4.307878(-f0) 




641 


2.327037(40) 


2.403445(40) 


1.208407(+0) 


2.087263(-l) 


642 


1.997378(40) 


-8.126205(-1) 






643 


1.079113(-1) 


1.762409(-2) 


1.096601(-2) 


3.058818(-3) 


b'44 


2.327409(40) 


6.901582(-1) 


-2.158431(-1) 


-1.084117(-1) 



640 = max (640, 1.0) 

041 = O41 

644 ~ b'44 
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a 




P 1 r] 


&46 


2.214315(- 


hO) 


-1.975747(+0) 


648 


5.072525(- 




1.146189(+1) 6.961724(+0) 1.316965(+0) 


649 


5.139740(- 







b45 = 1.0 - (2.47162p - 5.401682p2 + 3.247361p^) 
645 = 1-0 for p< 0.0 

6,. = -1.064a log (^) 

647 = 1.127733p + 0.2344416p2 - 0.3793726p3 





Q 


P 


7 


n 






1.125124(+0) 


1.306486(+0) 


3.622359(+0) 


2.601976(+0) 


3.031270(-1) 


652 


3.349489(-l) 


4.531269(-3) 


1.131793(-1) 


2.300156(-1) 


7.632745(-2) 




1.467794(+0) 


2.798142(+0) 


9.455580(+0) 


8.963904(+0) 


3.339719(+0) 


654 


4.658512(-1) 


2.597451(-1) 


9.048179(-1) 


7.394505(-l) 


1. 607092 (-1) 


655 


1.0422(+0) 


1.3156(-1) 


4.5000(-2) 






&56 


1.110866(+0) 


9.623856(-l) 


2.735487(+0) 


2.445602(+0) 


8.826352(-l) 


&57 


-1.584333(-1) 


-1.728865(-1) 


-4.461431(-1) 


-3.925259(-l) 


-1.276203(-1) 



651 - 0.1343798C'' 
6'53 + 0.4426929C'' 
min (0.99164 - 743.123Z2-*3 ^ ^^^.^ 

656 +0.1140142C^ 

657 - 0.01308728C^ 

Note that x(n) for some number x represents x x 10". 
A blank entry in a table implies a zero value. 



This paper has been produced using the Royal Astronomical Society /Blackwell Science KTgiX style file. 
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